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I.  INTRODUCTION 


Segmentation  techniques  are  among  the  most  important  considerations  in  the 
development  of  the  automated  image  processing  systems.  Two  related  algorithms  using 
2-D  linear  prediction  models  and  the  Karhunen-Loeve  transformation  for  multichannel 
and  color  image  segmentation  are  developed  and  compared  in  this  thesis. 

The  purpose  of  segmentation  is  to  partition  an  image  into  a  set  of  simpler 
homogeneous  regions.  The  regions  may  consist  of  different  gray  level,  different 
textures,  colors,  etc.  In  some  cases  an  '  image  '  may  consist  of  several  spectral 
components.  For  example,  a  color  image  consists  of  three  separate  signals  representing 
the  intensity  of  red,  green,  and  blue  as  a  function  of  spatial  position.  If  we  represent 
each  of  these  signals  by  functions  F^.  (n,m),  F,^  (n,m),  and  F^  (n,m),  the  image  is 
represented  by  a  vector  quantity 


i  (n,m)  - 


■F^(n,m) 
Fb  (n.m) 
.Fg(n,m)  , 


(1.1) 


where  n  and  m  represent  spatial  coordinates.  We  call  such  an  image,  consisting  of 
more  than  one  two  dimensional  (2*D)  signal,  a  multichannel  image. 

In  this  thesis,  a  method  based  upon  linear  prediction  is  evaluated  experimentally. 
This  method  has  been  developed  [Ref.  1}  for  monochrome  images  and  extended  to 
color  images  [Ref.  2].  That  method  uses  maximum  likelihood  {ML)  and  maximum  a 
posteriori  {MAP)  estimation  to  segment  multichannel  images  into  regions  of  similar 
textures.  The  linear  prediction  is  a  filtering  of  the  multichannel  image  to  estimate  the 
gray  level  at  a  particular  spatial  coordinate  from  the  gray  levels  at  neighboring 
positions.  It  is  implemented  as  a  2-D  linear  filtering  operation.  The  algorithm  uses  a 
previously-determined  set  of  parameters  corresponding  to  the  mean  of  the  data  in  each 
channel,  the  covariance  matrix  of  the  prediction  error,  and  the  weighting  coefficients  of 
the  estimation  filter  for  each  specific  texture  type. 

The  method  discussed  above  is  compared  to  a  variant  of  this  method  based  on 
the  Karhunen-Loeve  {K-L)  transformation.  The  K-L  transformation  allows  the  several 
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components  of  an  image  to  be  combined  into  a  single  image  that  retains  most  of  the 
energy  in  the  original  image.  Hunt  and  Kubler  [Ref.  3]  found  that  for  image 
restoration,  Karhunen-Loive  transformation  followed  by  single  channel  image 
processing  worked  nearly  as  well  as  multichannel  image  processing.  It  was  desired  to 
sec  if  the  Karhunen-Loive  transformation  would  be  equally  effective  for  segmentation. 
In  this  part  of  the  work,  the  K-L  transformation  has  been  developed  to  reduce  the 
3-channel  color  problem  to  a  1-channel  problem  and  the  segmentation  was  performed 
for  the  one-channel  image.  Karhunen-Loive  transfonnation  is  based  upon  the 
statistical  characteristics  of  an  image.  The  advantage  of  this  approach  is 
computational  savings;  only  about  one  ninth  the  number  of  computations  is  required 
for  this  method. 

The  remainder  of  this  thesis  is  organized  as  follows.  Chapter  II  discusses  the 
model  and  the  algorithms  used  to  perform  the  multichannel  image  segmentation 
employing  techniques  of  linear  prediction  [Ref.  4j.  A  general  class  of  linear  filtering 
models  for  texture  is  first  presented.  An  algorithm  is  then  developed  to  estimate  the 
filter  parameters  from  a  multichannel  image.  Then,  the  multichannel  image 
segmentation  algorithm  is  described. 

Chapter  III  presents  the  models  and  the  algorithms  to  perform  the  Karhunen- 
Lodve  transformation  and  one-channel  image  segmentation.  First,  the  algorithms  for 
determining  the  eigenvectors  and  the  eigenvalues  of  the  correlation  matrix  are 
developed.  Then,  the  transformation  using  a  3-channel  image  is  presented.  Finally, 
the  one-channel  image  segmentation  algorithm  is  discussed.  The  examples 
demonstrating  the  application  of  the  segmentation  methods  for  color  images  are 
presented  and  compared  in  Chapter  IV.  Chapter  V  has  the  conclusions  about  the 
multichannel  image  segmentation  and  one-channel  image  segmentation. 

In  Appendix  A,  the  Relaxation  Method  is  described  briefly  as  an  alternative  to 
the  maximum  a  posteriori  region  estimation  for  monochrome  images.  Results  are 
compared  with  the  MAP  method.  In  each  of  Appendices  B  and  C,  the  description  and 
use  of  the  computer  programs  for  the  multichannel  image  segmentation  and  the  one- 
channel  image  segmentation  are  presented.  The  computer  program  for  multichannel 
image  segmentation  is  contained  in  Appendix  D.  The  Karhunen-Loeve  transformation 
and  one-channel  image  segmentation  computer  programs  are  contained  in  Appendix  E. 
These  programs  are  written  in  FORTRAN,  compiled  using  Version  4.2  under  the  VAX 
/VMS  Version  4.1  operation  system. 
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n.  IMAGE  SEGMENTATION  MULTICHANNEL  FILTERING 


In  this  chapter,  a  multichannel  image  segmentation  algorithm  based  upon  a  2-D 
linear  filtering  model  is  presented.  The  multichannel  images  used  in  this  work  are  color 
images  with  three  channels  representing  the  red,  green,  and  blue  components.  The 
linear  filtering  model  is  used  to  develop  approximate  expressions  for  the  multivariate 
probability  density  functions  in  terms  of  the  filter  error  residuals  for  the  entire  set  of 
points  representing  an  image.  The  density  expressions  are  used  in  the  formulation  and 
solution  of  the  multichannel  image  segmentation  problem.  It  is  assumed  that  the 
multichannel  images  contain  multiple  regions  of  homogeneous  texture.  This  is  found  to 
be  the  case  in  dealing  with  aerial  photographs  of  natural  terrain. 

The  problem  of  multichannel  image  segmentation  is  addressed  as  an  estimation 
problem  for  two  regions  of  texture.  Maximum  likelihood  (ML)  and  maximum  a 
posteriori  {MAP)  region  estimates  using  a  Markov  random  field  to  model  region 
transitions  are  developed. 

This  chapter  consists  of  two  sections.  The  first  section  describes  the  linear 
filtering  model  and  develops  an  expression  for  the  image  probability  density  function  in 
terms  of  filter  error  residuals.  The  last  section  deals  with  the  algorithm  that  is 
developed  for  the  texture  estimation  and  the  multichannel  image  segmentation.  The 
results  of  texture  estimation  and  image  segmentation  are  presented  in  Chapter  IV. 

A.  MODEL  DEVELOPMENT 

In  this  section,  a  2-D  multichannel  autoregressive-moving  average  (ARMA)  model 
[Ref.  5]  is  first  discussed.  Then,  we  concentrate  on  the  multichannel  autoregressive 
(AR)  model  with  Gaussian  white  noise  inputs  [Ref  6].  The  development  parallels  that 
in  [Ref  1].  A  multichannel  image  is  represented  by  a  vector  signal  f*’  (n,m)  where 
(n,m)  are  spatial  coordinates  and  the  superscript  h  is  an  index  representing  the  texture 
type.  A  2-D  multichannel  image  is  shown  in  Figure  2. 1  .  A  texture  of  type  h  is  then 
modeled  in  general  by  a  multichannel  ARMA  process  defined  by 


(2.1) 


P-1  Q-1 

t  (n,m)  -  ^  (^.m) 

i-Oj-0  P 

(i,j)  *  (0,0) 


(n.m)  -  t  (n.m)  +  (2.2) 

for  h  *  0,1,  n  =  1,....,N,  m  =  1,....,M,  where  and  Bjj**  are  set  of  filter  weighting 
coefficient  matrices  of  size  K  by  K.  (n.m)  are  a  set  of  independent  identically 
distributed  zero-mean  random  variables,  £  is  a  constant  representing  the  mean  value  of 
the  image,  and  P  is  finite-extent  mask,  covering  the  filtered  points.  (n,m),  (n,m), 

and  are  vectors  of  size  K,  the  number  of  channels  in  the  image.  The  matrices  A*'.,  are 
the  key  parameters  in  the  linear  model.  For  a  first  quadrant  filter  with  a  P  x  Q  region 
of  support,  there  are  PQ  A^j  matrices  with  A*‘qq  =  I,  an  identity  matrix. 

For  the  auto-regressive  (AR)  or  all-pole  model  we  have  =  6^.. ’and  so  that 
Equation  2. 1  reduces  to 


P-1  Q-l 

t  (n.m)  -  -  II  A^ij  (n-i,m-i)  +  U  (n.m)  (2.3) 

i=0j  =  0 
(i.j)  *  (0,0) 

If  the  vectors  f,  w,  and  g  represent  an  ordered  set  of  the  corresponding  image  points, 
then  Equation  2.2  and  Equation  2.3  are  written  in  a  matrix  formulation  as 

A(f-  g)  -  W-Aofg  (2.4) 

where  A  and  A^  are  matrices  whose  nonzero  elements  are  derived  from  the  terms  A.,  in 

0  ij 

Equation  2.1  and  fg  represents  a  set  of  boundary  conditions  with  support  outside  of  the 
regions.  Since  the  terms  H  (n.m)  are  independent  with  probability  density  function 
(PDF)  p^,  One  can  solve  Equation  2.4  for  W  and  express  the  multivariate  probability 
density  function  for  the  image  conditioned  on  the  boundary  values  as 
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Figure  2.1  2-D  Multichannel  Image  Model. 


(2.5) 


Pf|  fo<f'^o)-  ^,Pw<A<f-*)  +  A^fo) 

•  n  (  £  (n,m)) 

(n,m)  e  R 

where  the  notation  E(n4n)  is  used  to  represent  the  ordered  components  of  the  vector 
A(  f-g)  Aq  Tq  .  If  the  boundary  conditions,  ,  are  temporarily  ignored,  then 


£  (n,m)  *  A  a  -  g) 


(2.6) 


and  the  terms  £  (n,m)  in  Equation  2.5  are  computed  from 


P-1  Q-l 

t  (n,m)  -  5:  S  A**.,  t  (n-i.m.j)  (2.7) 

i-Oj-0 
(i,j)  *  (0,0) 

The  filter  of  Equation  2.7  which  computes  £  **  (n,m)  from  £**  (n,m)  is  referred  to  as 
the  'prediction  error  filter'.  One  can  think  of  Equation  2.7  as  producing  an  estimate  or 
prediction  £  (n,m)  of  the  image  at  point  (n,m)  and  then  forming  an  error  £  (n,m)  as  a 
difference  £  (n,m)  -  ^  (n,m) .  This  process  is  known  as  2-D  linear  prediction  and  is 
fundamental  to  performing  multichannel  image  segmentation. 

B.  ALGORITHM  DEVELOPMENT 

In  the  multichannel  segmentation  problem,  it  is  assumed  that  the  image  consists 
of  multiple  connected  regions  of  known  texture  types,  but  that  the  region  boundaries 
and  the  number  of  regions  are  not  known.  The  segmentation  of  the  image  is  treated  as 
a  supervised  learning  problem,  since  the  regions  are  considered  to  consist  of  known 
texture  types.  In  this  section,  the  multichannel  image  segmentation  algorithm  for 
textured  images  is  discussed. 
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An  overview  of  the  method  is  as  follows.  Given  a  multichannel  image  of  each 
texture,  filter  parameters  are  estimated  by  computing  the  covariance  matrix  from  a  set 
of  data  and  solving  the  Normal  equations  corresponding  to  the  model  of  Equation  2.3  . 
In  this  case  the  'correlation  method  of  linear  prediction  is  used  to  compute  the 
covariance  matrix.  The  filter  parameters  are  derived  from  a  statistical  analysis  of  the 
textured  images,  because  the  image  model  discussed  in  the  previous  section  is  based 
upon  statistical  properties. 

Once  the  filter  parameters  are  known  the  filters  are  used  to  perform  the 
segmentation.  The  filter  weighting  coefficients  are  used  to  calculate  the  prediction 
errors  E**  (n,m)  of  two  textures  (n,m).  Then,  a  maximum  likelihood  (ML)  region 
estimate  is  developed  using  the  prediction  errors  and  the  covariance  matrices  for  image. 
The  ML  estimate  is  used  as  a  basis  to  determine  an  approximate  maximum  a  posteriori 
(MAP)  region  estimate.  The  MAP  region  estimation  utilizes  an  underlying  Markov 
structure  for  the  region  statistics  to  produce  an  accurate  segmentation. 

1.  Filter  Parameter  Estimation  Method 

The  prediction  error  filter  is  a  finite-extent  impulse  response  (FIR)  or 
nonrecursive  filter  with  selectable  mask  size  and  quarterplane  region  of  support.  It  is 
always  stable.  The  inverse  of  the  filter  used  in  Equation  2.3  is  an  alLpole  filter  (the  AR 
filter).  The  filter  parameter  estimation  problem  requires  calculating  ,  'LJ',  and  A** 
by  statistical  analysis  of  data  in  an  estimation  window  containing  the.  desired  texture, 
where  the  quantity  if**  is  a  njean  vector  of  the  average  gray  level  of  the  image  in  each 
of  the  channels,  and  is  a  covariance  matrix  of  the  multichannel  white  noise 


2:^  -  E  ( (n,m)  (^  (n,m))T  ] 


(2.8) 


the  term  E*’^  is  also  refered  to  as  the  prediction  error  covariance  matrix,  since  in  a 
linear  prediction  problem  it  represents  the  covariance  of  the  quantities  g  (n,m)  defined 
in  Equation  2.7  .  Since  E**^  is  not  in  general  a  diagonal  matrix,  we  see  that  the  white 
noise'  is  uncorrelated  within  each  channel  but  correlated  between  channels. 

The  covariance  of  the  multichannel  white  noise  (the  prediction  error 
covariance)  and  the  filter  coefficients  will  be  obtained  by  estimating  the  correlation 
function  of  the  image  and  solving  a  set  of  Normal  equations  as  discussed  below.  The 
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correlation  function  itself  is  estimated  from  data  in  a  window  containing  a  sample  of 
the  desired  texture.  Two  estimation  windows  are  depicted  in  Figure  2.2  for  two 
different  textures  in  the  multichannel  image.  The  reader  should  realize  that  the 
estimation  windows  do  not  have  to  come  from  the  image  to  be  segmented;  they  can  be 
selected  from  any  image  containing  the  same  type  of  texture. 


a.  A/e<in  Vector  Estimation 

In  order  to  model  the  multichannel  image  by  Equations  2.2  and  2.3  the 
mean  vector  of  the  multichannel  image  has  to  be  estimated.  Knowing  ^  and  selecting 
the  stationary  estimation  windows  of  two  desired  textures  of  Zjn>  mean  2-D 

multichannel  image,  £,  appearing  in  Equation  2.3  can  be  obtained  by  subtracting  the 
mean  vectors  from  the  multichannel  image.  Thus,  Equation  2.2  becomes 


t  (n,m)  -  (n,m)  -  it 


(2.9) 
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wlicrc  corresponds  to  g*'  in  liquation  2.2,  and  the  term  fjf'  in  Equation  2.9  arc 
computed  Prom 

it  '  M''2  (2.10) 

.^*•3  . 

where 

1  X2  Y2 

- ^ —  I S  <2. 1 1) 

n' m'  n=X|m-Yi 

is  the  mean  estimate  for  the  k*  channel  of  the  h*  texture  image.  The  limits  Xj  , 
Xj  ,  Yj  ,  Yj  represents  the  edges  of  the  window  which  is  of  size  N'  by  M'.  Therefore 
n'  =  X2  -  Xj  +  1  ,  m'  -  Yj  -  Y,  +  1,  and  0  ^  Xj  <  Xj  ^  N,  0  ^  Yj  <  Y^ 
^  M  .and  h  represents  the  two  textures  of  the  multichannel  image.  All  variables  used 
in  Equation  2.11  are  depicted  in  Figure  2.3  . 


Figure  2.3  Selection  of  the  Size  of  the  Windows  for  Mean  Vector  Estimation. 
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b.  Correlation  Function  Estimation 

The  correlation  function  of  the  zero  mean,  2-D  multichaimel  signal  has  to 
be  calculated  in  order  to  estimate  the  multichannel  white  noise  covariance  or  prediction 
error  covariance,  and  the  filter  weighting  coefficients.  A**..  .  The  theoretical  2-D 

W  IJ 

matrix  correlation  function  for  lag  "  i,j '  is  given  by 

R**  (i.j)  -  (R**  (-i.  (2.12) 

-  E  [£*>  (n,m) .  (  (n-i,m-j))’^  ] 

and  can  be  estimated  from  the  multichaimel  signal  by 
1  n2  m2 

R**  (i.i)  - - 1 1  (n,m)  (  (n-i,  m-j)^  (2. 1 3) 

*  »  » 

N  M  n*n|  m=mj 

where  R**  (i,j)  is  a  matrix  of  size  K  by  K,  and  nj  ,  nj  ,  mj  ,  nij  are  defined  by 

0  ^  nj  *'  max(Xj  ,X|  +  i )  <  n2  =  minCXj  ,  Xj  +  i )  i  N,  and 
0  ^  mj  “  max(Yj  ,  Yj  +  j)  <  m^  “  min(Y2  ,  Yj  j )  i  M. 

This  matrix  correlation  function  is  used  to  form  a  larger  block  Toeplitz  covariance 
matrix  which  is  used  to  estimate  the  filter  parameters.  This  is  discussed  next. 

c.  Filter  Coejfjficients  .and  Prediction  Error  Covariance 

The  prediction  error  filter  weighting  coefficients  and  the  prediction  error 
covariance  must  satisfy  a  set  of  linear  equations  known  as  the  Normal  Equations  when 
the  multichannel  image  is  represented  by  the  model  in  Equation  2.3  . 

Normal  equations  corresponding  to  Equation  2.3  can  be  written  as 

[R]  .  [A]  -  IS]  (2.15) 

where  R  is  the  correlation  matrix  for  the  signal,  A  is  an  appropriately  ordered  matrix 
of  the  filter  coefficients,  and  S  contains  a  single  non-zero  block  which  is  the 
prediction  error  covariance.  The  matrix  R  has  three  levels  of  partitioning  and  for  any 
rectangular  region  of  support  is  block  Teoplitz  with  block  Toeplitz  blocks. 

For  a  first  quadrant  filter  with  a  P  x  Q  region  of  support,  The  Normal 
equations  that  define  Equation  2.15  have  the  specific  form 
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(2.16) 


■R(0)  R(.l)  .  .  .R(.p+1)' 
R(l)  R(0)  .  .  .R(.P+2) 

r  ao  *1 

A' 

.gO 

0 

.R(P-1)  R(P-2).  .  .  R(0)  . 

AP-'. 

0 

m 

where 


R(i) 


R(i.O)  R(i,-1)  . 
R(i,l)  R(i.O)  . 


•  R(4-Q+  1) 

.  R(i.-0+2) 


LR(i,Q-l)  R(i,Q-2)  .  .  R(i,0) 


(2.17) 


and  where  R(i,j)  is  the  matrix  correlation  function  described  in  Equations  2.12  and  2:13 
The  quantities  A‘  and  are  deflned  by 


A*  - 


(2.18) 


A: 


.Q-i  J 


and 


S  - 


(2  20) 


L  0 
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with  Aqq*  I  and  where  A-  and  the  partitions  of  S®  are  matrices  of  order  K,  the 
number  of  channels  of  image. 

2.  Multichannel  Segmentation  Method 

An  overview  of  the  segmentation  is  as  follows.  By  using  a  Gaussian 
probability  density  for  the  white  noise  in  Equation  2.5  ,  one  can  develop  explicit 
estimates  for  the  density  functions  in  terms  of  the  prediction  errors  S  (n.m) .  From  this 
one  can  form  the  conditions  for  ML  and  MAP  estimation  of  the  regions  in  the  image. 
The  theory  leading  to  the  estimates  is  explained  below.  A  brief  intuitive  explanation  of 
the  process  is  given  here. 

As  mentioned  earlier  the  prediction  error  filters  in  Equation  2.3  can  be 
considered  as  predicting  the  intensity  of  a  pixel  in  each  channel  from  data  in  the  region  , 
adjacent  to  the  vector  of  pixels.  The  prediction  errors  (n,m)  are  the  outputs  of  the 
filters.  In  the  segmentation  algorithm  the  prediction  error  is  normalized  in  an 
expression  involving  the  corresponding  prediction  error  covariance.  These  normalized 
errors  are  compared  in  an  appropriate  formula  to  obtain  the  ML  region  estimate. 
When  an  area  of  texture  is  processed  by  a  filter  that  is  not  matched  to  the  texture,  the 
normalized  prediction  error  can  be  expected  to  be  high.  When  the  same  area  is 
processed  by  the  filter  that  is  matched  to  the  texture,  the  prediction  error  can  be 
expected  to  be  low. 

The  ‘maximum  likelihood  region  estimate,  ML(n,m),  of  texture  class  is 
achieved  based  on  the  prediction  enor  covariance  and  the  error  estimates  of  the  two 
desired  textures  for  each  pixel  in  the  multichannel  image.  The  .ML(n,m)  region 
estimate  assigns  pixels  to  texture  types  without  regard  to  the  assignments  of  the 
adjacent  pixels.  Then,  the  ‘maximum  a  posteriori  region  estimate,  MAPfn.m),  of 
texture  class  is  achieved  for  a  pixel  and  a  desired  number  of  adjacent  pixels  of  two 
textures  of  the  ML  region  estimation  result.  The  MAP  region  estimation  uses  the 
Markov  model  that  refers  to  the  above  description.  The  form  of  the  Markov  model 
and  ML  and  MAP  region  estimates  are  presented  in  detail  below. 

a.  Maximum  Likelihood  Region  Estimation 

It  is  supposed  that  a  multichannel  image  has  many  regions,  but  that  each 
region  contains  only  one  or  another  of  two  texture  types.  Given  these  regions,  one  can 
write  the  Equation  2.5  as 

P(f  I  R|)  -  n  p^  iS  (n,m))  i  -  1 . q  (2.21) 

hi 

(n.m) 
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where  the  is  the  probability  density  function  for  the  white  noise  source  of  type  hj 
within  region  Rj  .and  q  is  the  number  of  regions.  When  the  white  noise  term  is 
Gaussian  with  density  function  (mean  0  and  covariance  ) 


\ 


( 


NM/2 
2it)  1 


1.1 


fli 


expi  - 


(2-22) 


then,  taking  minus  twice  the  log  of  Equation  2.21  and  applying  Equation  2.22  ,  we 
obtain  for  an  N  by  M  pixel  multicharmel  image 


-  21np(f  (Rj.Rj . \) 

-I  (n.m)lT  [Z\  r  '  \M\  (n.m)l  +  ln|L;i  +  .. 

(n,m)  c  R, 

+  I  (Iff’h  r  ^  (n.m)]  +  ln|L%  |) 

(n,m)  e  R^  —  NM  ln2Jt 

(2.23) 

q 

-II  (tfh  r  '  (ffh  (n.m)]  +  lnir^  I) 

i*l(n,m)6R. 

-  NM  ln2Jt 


For  maximum  likelihood  estimation,  the  number  of  regions  q  and  the  regions 
themselves  are  considered  to  be  deterministic  parameters  of  the  density  function.  An 
ML  estimate  for  these  parameters  is  obtained  by  choosing  values  that  maximize 
Equation  2.21  or,  minimize  Equation  2.23  .  Since  NMln2R  is  constant  value,  the 
Eqtiation  2.23  is  minimized  if  every  point  (n,m)  in  the  multichannel  image  to  a  region 
Rj  of  type  hj  such  that  the  term  in  brackets  is  minimum.  Thus,  one  can  write  a  ML 
region  estimation  for  two  textures  as 


ML 


0 

(n,m)  > 
1 


MLo(n,m) 


(2.24) 


where 


4 
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MI^  (n.m)  -  (n.m)lT  J*’  [£*>  (n.m)l+lnlL»»^l 


(2.25) 


for  h  ■  0  ,  1,  where  the  number  above  or  below  the  inequality  indicates  the  region 
class  to  which  the  pixel  (n,m)  is  assigned.  When  the  class  is  1,  the  ML(n,m)  is  assigned 
I  for  a  white  pixel.  Otherwise  ML(n,m)  is  assigned  0  for  a  black  pixel.  Since  the  ML 
region  estimate  assigns  pixels  to  black  and  white  without  regard  to  the  assignments  of 
adjacent  pixels,  this  algorithm  produces  a  number  of  false  assigiunencs  and  a  somewhat 
spotty"  result. 

b.  Maximum  A  Posteriori  Region  Estimation 

The  'maximum  a  posteriori'  (MAP)  region  estimation  utilizes  the  Markov 
model  to  describe  the  occurance  of  regions  in  the  image.  The  combination  of  the  linear 
filtering  model  with  the  Markov  model  results  in  an  algorithm  to  achieve  a  MAP 
region  estimation.  For  MAP  estimation  the  regions  are  considered  to  be  random 
quantities,  and  we  maximize  the  probability  for  a  given  set  of  regions  conditioned  on 
our  observation  of  the  multichannel  image.  From  Bayes  rule,  the  a  posteriori 
probability  can  be  written  as 


Pr[K^  ,Rj . \\l] 


P(f  . . )  *  ^^^1  »  ^2  * . *^Q  ^ 

p(f) 


(2.26) 


Since  the  denominator  of  Equation  2.26  does  not  depend  on  the  regions,  maximum  a 
posteriori  (MAP)  estimate  for  the  regions  can  be  obtained  if  the  R.  are  chosen  to 
maximize  the  numerator 


piE  IRi  .Rj . \  ) .  Pr  (Ri  .Rj . \  1  (2.27) 

One  can  define  the  ”  state  "  s(n,m)  of  point  (n,m)  as  the  region  type  to 
which  that  point  has  been  assigned.  In  our  development,  the  number  of  region  types  is 
assumed  to  be  2.  Since  the  set  of  all  possible  state  assignments  for  points  in  the  image 
is  one-tO'One  with  the  set  of  all  possible  divisions  of  the  multichannel  image  into 
regions,  the  region  estimation  problem  can  be  viewed  as  one  of  estimating  the  states  of 
the  points.  It  is  assumed  that  the  state  of  a  point  is  stochastically  dependent  on  some 
adjacent  set  of  states  ^  in  a  s>mmetric  support  region,  as  shown  in  Figure  2.4  . 
Since  S  represents  a  chosen  set  of  state  assignments  for  all  points  in  the  multichannel 
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image,  Pr  [  S  ]  denotes  the  joint  probability  that  the  points  in  the  image  take  on  a 
chosen  set  of  state  assignments.  The  support  set  defines  a  neighborhood  structure  i.e. 
all  elements  in  the  support  set  are  neighbors  of  each  other.  It  can  be  shown  that  if  the 
set  of  states  S  is  a  Markov  random  field,  then  the  probability  of  S  can  be  factored  as  a 
product  of  terms  of  functions  depending  only  on  the  "cliques"  of  the  support  set  ^  . 
The  cliques  are  defined  as  groups  of  points  such  that  each  set  of  points  are  neighbors 
of  each  other  according  to  the  support  set.  For  a  Markov  random  field,  the  probability 
of  S  can  be  written  as 


o  o  o 

U  V  W 

e  o  e 

t  S  t' 

o  a  o 

w'  v'  u' 

I _ ^ _ 

^n,  m 

Figure  2.4  State  support  region  of  the  point  s  for  MAP  estimation. 

Pr  (  S  ]  -  n  Pr  (  s(n,m)  I  S„  „  ]  (2.28) 

‘  ^  (n,m)  '  '  ’  n,mi 

where  the  terms  in  the  product  are  additive  functions  defined  on  the  cliques  of 
and  the  product  is  over  all  cliques  in  S.  One  simple  acceptable  form  for  the  terms  in 
the  product  is 


Pr  [  s(n,m)  I  ^  exp  ( s(n,m)  {o  +  Pj  (t  +  f)  -t-  p^  (v-»- v  ) 

+  Yj  (u+u')  +  Y2  (w+w  )}]  (2.29) 

where  t  -  s(n-l,m),  t'  -  s(  n+ 1,  m),  ^  is  the  set  of  states  as  snown  in  Figure  2.4  , 

and  D  is  a  normalizing  constant.  One  particular  selection  of  the  parameters,  namely  a 
-  -  4,  Pj  -  Pj  *  Yj  *  Y2  *  *0  2  particularly  simple  algorithm.  In  this 

case,  we  have 
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(2.30) 


-  exp  (  (  s(  i,  j )  -  1  /  2  )) 

D  (i.j)  e  S„.  „ 


Pr  [  0  I  S„  „  ]  -  - 

'  n,  m  J 


(2.31) 


The  second  term  in  the  numerator  of  Equation  2.26  can  be  replaced  by 
Pr(S] .  Thus  maximizing  Equation  2.26  is  equivalent  to  minimizing 


-  21n  p(  f  I  R, . R  )  -  21n  Pr  [  S  ] 


(2.32) 


One  can  define  the  MAP  region  estimate  by  combining  Equations  2.22  , 
2.26 . 2.28  .  and  2.32  as 


^jS'  (n.m)) 


+  ln|i:V|-2tal>r|l|S„  „,!< 

1 


(n,m)l*  1  [f  (n.m)l 


(2.33) 


+  ln|i:Vl-21nPr[0|S„^^] 


For  Pr  I  h  I  Sjj  ijj  ]  in  the  form  of  Equation  2.30  and  Equation  2.31  computing  the 
terms  -  2  In  Pr  (hlS^j^  ]  is  equivalent  to  counting  the  number  of  pixels  in  that 
have  value  'h'  and  dividing  by  the  total  number  of  pixels  in  ^  ,  and  multiplying  by 
an  appropriate  normalizing  factor,  KS.*  A  larger  state  support  region  is  depicted 
in  Figure  2.5  as  an  example.  The  side  of  the  must  be  an  odd  number  .  Although 
it  does  not  necessarily  lead  to  the  true  MAP  estimate  we  find  it  convenient  in  practice 
to  maximize  Equation  2.33  term  by  term.  That  is,  we  require  to  use  Equation  2.33 
without  sum.  Equation  2.33  can  be  solved  by  iteration  using  the  maximum  likelihood 

1  The  normalizing  factor  results  because  the  quantities  a  ,  B,  ,  B,  ,  Yi  ,and  y,  in 
mnction'  2.29  can  be  scaled  arbitrarily  and  stilf  result  in  a  liguimate  *  probabftity 
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state  estimates  obtained  from  liquation  2.24  as  an  initial  set  of  states  for  MAP  region 
estimate. 
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Figure  2.5  A  Set  of  States  1  's,  O's  adjacent  to  S  for  MAP  Region  Estimate. 


1  he  computational  requirements  of  the  MAP  region  estimation  are  reduced  by  storing 
the  differences, 


MLD(n,m)  -  ML|  (n,m)  -  MLQ(n.m)  (2.34) 

incurred  during  the  calculation  of  the  ML  region  estimate.  Substituting  Equation  2.34 
in  Equation  2.33  gives 

MLD(n,m)  “  2  In  Pr  [  1  |  S„  )  +  2  In  Pr  [  0  1  S„  „,  ]  <  0  (2.35) 

0 

At  each  iteration  terms  Fr  (  h  I  j  are  evaluated  based  on  the  values  of  the  states 
at  the  previous  iteration.  For  our  particular  method  of  selecting  Pr  (  h  |  S^^  ^  ]  , 
Equation  2.35  can  be  expressed  as 


MLD(n,m) 


-  KS 


2  (number  of  state  I  pixels)  -  (  number  of  pixels  in  S  )  +  1 

^  Df  fll  ^ 


number  of  pixels  in  S_ 

^  n  ,  m 


1 

<  0 
> 

0 


« 


(2.36) 
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There  are  two  importaitt  points  in  Equation  2.36  in  order  to  perform  the  maximum  a 
posteriori  image  segmentation  accurately.  The  value  of  the  convergence  factor,  KS, 
must  be  assigned  properly.  If  KS  is  assigned  too  small,  the  segmentation  may  not 
remove  improperly  classified  pixels.  On  the  other  hand,  if  KS  is  assigned  too  large, 
correctly  classified  pixels  could  be  changed.  In  addition,  the  size  of  must  be  large 
enough.  Otherwise,  false  assignments  produced  by  the  initial  ML  segmentation  may 
not  be  removed. 


24 


III.  ONE-CHANNEL 


In  this  chapter,  the  models  and  relevant  algorithms  are  presented  to  perform  the 
Karhunen-Lolve  (K*L)  transformation  [Ref.  3]  and  one-channel  image  segmentation 
utilizing  the  techniques  of  linear  prediction  [Ref  IJ.  Since  linear  prediction  techniques 
are  presented  in  detail  for  multichannel  image  segmentation  in  the  previous  chapter, 
this  discussion  concentratet  on  the  K-L  transformation.  The  K-L  model  and  algorithm 
are  first  developed  to  reduce  the  three-channel  color  problem  to  a  one-channel 
problem.  Then,  a  one-channel  segmentation  procedure  is  presented  that  is  based  on  the 
'same  model  previously  discussed.  The  results  of  the  K-L  transformed  one-channel 
image  segmentation  are  pres;-nted  and  compared  with  the  multichannel  image 
segmentation  in  Chapter  IV. 

A.  MODEL  DEVELOPMENT 

The  K-L  transformation  developed  in  this  section  is  based  on  the  statistical 
properties  of  an  image.  This  transformation  provides  an  energy  compaction  between 
channels  of  a  color  image.  That  is,  most  of  the  color  image  energy  is  compacted  into 
one  channel,  and  the  transformed  image  channels  are  uncorrelated.  If  the  multichannel 
image  and  transformed  multichannel  image  are  expressed  in  vector  form.  The  K-L 
transformation  is  given  by  [Ref  7) 

fi(n,m)  =[A]  f(n,m)  (3.1) 

where  £  (  n  ,  m  )  is  the  original  multichannel  image,  g  (  n  ,  m  )  is  the  transformed 
multichannel  image,  and  A  is  the  K-L  transformation  matrix,  whose  rows  are 
eigenvectors  of  the  between-channel  correlation  matrix  R  defined  by 

R  =  E[£(n.m)f’^(n.m)].  (3.2) 

The  betw'een-channel  correlation  matrix  of  the  transformed  image  is 

A“E[g(n,m)g^(n.m)]  (3.3) 

-  [A](R][Ar* 
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where  the  matrix  A  is  a  diagonal  matrix  of  the  form 


A  * 


•  X.J  0  0 

0  ^2  0 

.0  0  K3 


(3.4) 


and  the  represent  eigenvalues  of  the  between-channel  correlation  matrix.  The 
eigenvalues  are  ordered  such  that 


Xj  ^  Xj  ^  X3  ^  0 


(3.5) 


The  importance  of  this  property  is  that  each  eigenvalue  X^^  is  equal  to  the  variance  of 
the  k  *  channel  of  the  transformed  multichannel  image  whose  channels  are 
uncorrelated.  Then,  it  is  a  well  known  property  of  multivariate  statistics  that  the  total 
variability  of  the  color  image  has  the  form  [Ref.  3] 

3 

.  (3-^) 

k=  I 


which  relates  the  total  variability  to  the  decorrelated  component  variations,  Xj^  .  One 
can  observe  that  often  the  X^  values  have  a  wide  range  of  magnitudes,  and  the  first 
component,  Xj  ,  will  be  sufficient  to  approximate  Xj  with  only  a  small  percentage  of 
error.  This  becomes  the  key  idea  for  the  use  of  the  K-L  transformation.  Table  1  shows 
the  energy  distribution  between  the  transformed  color  image  channels  of  two  test 
images. 

Indeed,  in  a  3  channel  image  one  can  often  find  that  the  first  channel  of  the  K-L 
transformed  image  is  sufficient  to  account  for  99  percent  of  all  the  variability.  That  is, 
it  is  typical  to  find  that 

Xj  =  0.99  Xj  (3.7) 

The  Karhunen-Loive  transformation  provides  the  best  energy  compaction  [Ref.  8]  and 
the  advantage  of  this  transformation  is  computational  savings.  We  will  see  later  that 
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TABLE  1 

ENERGY  DISTRIBUTION  BETWEEN  TRANSFORMED  IMAGE 

ClIANNE^LS 


Percentage  of  Energy  in  Channels 

Qi 

02 _ 03 _ 

Image  1 

99.13 

Image  2 

99.15 

one-channel  image  segmentation*  achieved  by  processing  only  the  first  channel  of  the 
K-L  transformed  color  image  will  be  very  close  to  the  result  obtained  on  the  original 
image  with  a  multichannel  algorithm,  but  will  be  obtained  at  only  one  ninth  of  the 
computational  cost. 

B.  ALGORITHM  DEVELOPMENT 

In  this  section,  the  algorithm  to  perform  K-L  transformation  of  color  images  is 
presented.  Since  the  K-L  transformation  is  based  upon  the  between-channel  correlation 
matrix,  R  ,  of  the  color  image  ,  the  between-channel  correlation  matrix  is  first 
determined.  Then  the  eigenvectors  of  the  between-channel  correlation  matrix  are 
computed.  Finally,  the  K-L  transformation  is  formulated  using  the  transpose  of  the 
eigenvector  matrix. 

1.  Correlation  Function  Esthnation 

In  order  to  determine  the  K-L  transformation  of  Equation  3.1  ,  the  between- 
chatmel  correlation  matrix  of  the  color  image  must  first  be  estimated.  Our  estimate  for 
the  between-channel  correlation  matrix  is  given  by 


1  N-1  N-I 

R  - - 

N^  n-0  m-0 


(3.8) 
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where  the  image  size  is  N  by  N,  and  the  between-channel  correlation  matrix  size  is  K 
by  K. 

2.  Karhunen-Loeve  Transformation  Matrix 

Since  the  rows  of  the  K-L  transformation  matrix  are  the  eigenvectors,  E,  of 
the  between-channel  correlation  matrix,  the  eigenvectors  must  be  calculated  from  the 
between-channel  correlation  matrix  of  the  color  image.  The  K-L  transformation 
matrix  is  then  obtained  using  the  transpose  of  the  eigenvector  matrix.  That  is 


(3.9) 


« 

3.  Karhunen-Loeve  Transformation 

Since  the  K-L  transformation  matrix  satisfies  Equation  3.3  ,  then  the  K-L 
transformation  is  represented  by  Equation  3.1  with  A  is  given  by  Equation  3.9,  where 
eT  ,  i  -  1  ,  2  ,  3  are  the  eigenvectors  of  R.  More  explicitly,  the  components  of  S  are 
determined  from  the  components  of  £  at  any  pixel,(  n  ,  m  )  as 


■Qj  (  n  ,  m )  ■ 

F|  ( n,  m) 

a 

O' 

= 

F2  (  n ,  m  ) 

Qj  (  n  ,  m ) 

^ _ ^T  _ 

F3  (  n  ,  m ) 

(3.10) 


C.  ONE-CHANNEL  IMAGE  SEGMENTATION 

The  one-channel  segmentation  procedures  presented  in  this  section  are  based  on 
the  same  model  that  was  described  in  detail  for  the  multichannel  image  in  Chapter  II. 
A  single  channel  image  is  used  instead  of  a  three  channel  image.  That  is,  K  is  always 
equal  to  1  in  the  Equations  of  Chapter  II.  As  a  result  the  vector  and  matrix  quantities 
become  scalars  and  the  equations  are  simplified. 

In  order  to  model  the  single  image  by  Equations  2.2  and  2.3,  the  mean  of  the 
image  is  estimated  using  the  Equation  2.11  .  Then  the  correlation  function  is  estimated 
in  the  same  manner  as  in  Section  II-B.b  where  R(i,j)  is  a  scalar  instead  of  a  matrix. 
This  procedure  is  followed  by  estimation  of  the  filter  coefficients  and  the  prediction 
error  covariance.  The  equations  in  Section  II-B.c  are  then  used  to  determine  the  filter 
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coefficients,  A..  ,  and  the  the  prediction  error  covariance,  T  ,  now  a  scalar  value. 
Finally,  the  one-channel  segmentation  algorithm  is  applied.  The  method  is  the  same  as 
that  described  in  detail  for  the  color  images  in  Section  I1-B.2.  However,  instead  of 
Equation  2.24  and  2.33  the  following  simplified  relations  are  used  for  the  maximum 
likelihood  and  the  maximum  a  posteriori  region  estimates.  A  maximum  likelihood 
region  estimate  for  a  single  channel  image  of  the  pixel  is  given  by  [Ref  1] 


(  (n,m))2 


(n,m))'^  >  (  E“  (n.m))'^  ~ 

+  1”  (  E',  )  < -!— +  In  (  L”,  ) 


(3.11) 


and  the  maximum  a  posteriori  region  estimate  is  given  by 
(  E^  (n,m))2 

'-FT— +  ln(5:‘w)- 

0  **  w 

>  (E**  (n,m))2  „ 

to  (  S".  )  -  2  In  Pr  1  0  I  S„  1 


(3.12) 


where  E^  and  E°  are  the  result  of  applying  the  linear  predictive  filters  to  the  first 
channel  of  the  K-L  test  image. 
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IV.  RESULTS  AND  COMPARISON  OF  THE  METHODS 

In  this  chapter,  the  results  of  the  multichannel  image  segmentation,  the 
Karhunen-Locve  transformed  one-channel  image  segmentation,  and  the  K-L 
transformed  3-channel  image  segmentation  are  presented  and  compared. 

The  digitized  image  size  used  in  this  work,  is  128  by  128  pixels  with  gray  levels 
represented  on  a  scale  of  0  to  255  (8  bits).  A  digitized  color  photograph  of  a  rural  area 
containing  trees  (the  green  region)  and  fields  (the  yellow  region)  are  shown  in  Figure 
4.1  .  A  quarter-plane  filter  for  each  texture  class  (2^2  pixels)  was  designed  and 
applied  to  the  color  image  to  achieve  the  multichannel  image  segmentation.  A  state 
support  region  of  7  by  7  pixels  was  used  for  MAP  region  estimation.  The  results  of  the 
maximum  likelihood  (ML)  and  the  maximum  a  posteriori  (MAP)  segmentations  of 
Figure  4.1  image  are  shown  in  Figures  4.2  and  4.3  respectively.  The  segmentation 
results  show  the  field  regions  as  black  and  tree  regions  coded  as  white.  The  ML  result 
(Fig  4.2)  is  spotty,  but  the  true  tree  and  field  regions  are  distinguishable.  The  MAP 
segmentation  result  (Fig  4.3)  of  Figure  4.1  image  is  quite  clear.  The  MAP 
segmentation  was  not  able  to  remove  a  few  improperly  classified  points  in  the  left  side 
of  the  field  /egion.  Since  there  is  an  inherent  ambiguity  in  the  estimation  of  the  region 
boundaries  due  to  the  finite  size  of  the  masks,  the  edges  are  expectedly  somewhat 
.rough.  Figure  4.4  shows  another  color  image  of  a  rural  area  containing  trees  and 
fields.  Figures  4.5  and  4.6  show  the  results  of  the  ML  and  the  MAP  segmentation  of 
the  Figure  4.4  image.  The  ML  estimation  was  achieved  using  the  filters  designed  for 
the  image  of  Figure  4.1  .  The  result  of  the  ML  segmentation  (Fig  4.5)  is  again  quite 
spotty,  i.e.  Although  two  regions  are  discernible.  The  MAP  algorithms  segmented  the 
regions  quite  accurately  as  shown  in  Figure  4.6  .  The  MAP  estimates  presented  above 
converged  after  10  iterations  and  KS  was  assigned  to  100. 

In  the  results  presented  above  the  ML  procedure  produced  a  poor  result  with  a 
lot  of  false  regions.  This  is  due  to  the  lack  of  prior  information  about  region 
connectivity.  On  the  other  hand,  MAP  estimation  using  the  Markov  model  to  represent 
region  transitions  produced  results  that  was  quite  accurate. 

Figure  4.7  shows  the  first  channel  of  the  Karhunen-Loeve  (K-L)  transformed 
image  of  the  Figure  4.1  .  The  image  size  is  128  by  128  pixels  with  the  scaled  gray  levels 
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within  the  intensity  range  of  the  display.  The  result  of  the  one-channel  ML 
segmentation  of  Figure  4.7  is  depicted  in  Figure  4.8  .  The  one-channel  ML 
segmentation  result  is  quite  spotty,  but  the  two  regions  are  perceptually  discernible.  On 
the  other  hand,  the  MAP  segmentation  result  (Fig  4.9)  of  Figure  4.7  is  quite  smooth, 
although  there  are  a  few  incorrectly  classified  points  in  both  regions. 

The  first  channel  of  the  Karhunen-Loive  transformed  image  of  Figure  4.4  is 
shown  in  Figure  4.10  ,  and  Figures  4.11  and  4,12  show  the  results  of  the  one-channel 
ML  and  MAP  segmentations  of  Figure  4,10  respectively.  The  one-channel  ML 
segmentation  result  (Fig  4.11)  has  a  lot  of  spots,  but  both  segmented  regions  are 
distinguishable.  The  maximum  a  posteriori  segmentation  (Fig  4.12)  of  Figure  4.10  is 
quite  smooth,  but  again  there  are  a  few  incorrectly  classified  sets  of  pixels  in  both 
regions. 

The  results  of  multichannel  image  segmentation  and  the  Karhunen-Loeve 
transformed  one-channel  image  segmentations  were  presented  consecutively.  These 
results  show  that  the  ML  estimation  of  the  Karhunen-Loeve  transformed  single 
channel  image  is  much  more  spotty  than  the  ML  estimation  of  the  multichannel  image. 
But,  the  MAP  estimation  results  of  the  K-L  transformed  single  channel  images  are  as 

clear  as  the  MAP  segmentation  results  of  the  multichannel  images. 

* 

In  summary  the  results  presented  in  this  chapter  show  that  the  best  segmentation 
result  is  provided  by  the  multichannel  image  segmentation  method.  However  the 
seginentation  results  of  the  Karhunen-Loeve  transformed  single  charmel  images  are 
very  close  to  multichaimel  image  segmentation  results,  i.e.  Because  most  of  the  color 
image  energy  (99  percent)  is  compacted  into  the  first  ohaiuiel. 


Figure  4.1  Two-Texture  Color  Image. 
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Figure  4.5  ML  Region  Estimation  of  Figure  4.4  Image. 


Figure  4.6  MAP  Region  Estimation  of  Figure  4.4  Image. 


V.  CONCLUSIONS 


The  segmentation  of  terrain  images  is  an  important  part  of  image  analysis 
methods  for  military  and  civilian  applications.  The  work  in  this  thesis  utilized  a  2-D 
stochastic  linear  filtering  model  and  compared  algorithms  for  multichannel  image 
segmentation  of  color  images.  Two  levels  of  structure  were  used  for  multichannel 
segmentation  development.  The  fundamental  structure  based  on  the  linear  filtering 
concepts  represents  the  texture  in  local  regions  of  terrain.  Superimposed  on  this 
structure  is  a  Markov  random  field  that  describes  transitions  from  one  region  type  to 
another.  The  segmentation  was  considered  as  a  region  estimation  problem  and 
maximum  likelihood  and  maximum  a  posteriori  region  estimatation  methods  were 
developed.  The  ML  region  estimation  produced  a  spotty  result,  but  the  MAP  region 
estimation  produced  quite  accurate  results  for  the  multichannel  and  single  channel 
image. 

The  other  piece  of  work  developed  in  this  thesis  was  the  Karhunen-Loeve 
transformation  model  that  based  on  the  statistical  characteristics  of  color  image.The 
one-channel  image  segmentation  was  then  applied  to  the  first  channel  of  the 
Karhunen-Loeve  transformed  color  image  to  sec  the  effectiveness  of  the  K-L 
transformation  for  segmentation. 

We  observed  that  multichannel  image  segmentation  results  were  quite  accurate. 
Similarly  the  results  of  the  K-L  transformed  one-channel  image  segmentation  were  very 
smooth.  In  summary  the  results  of  both  segmentation  methods  were  very  close  to  each 
other,  and  the  K-L  transformation  is  very  effective  for  segmentation. 
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APPENDIX  A 
RELAXATION  METHOD 


The  Relaxation  Method  algorithm  utilizes  a  set  of  iterative  numerical  techniques 
to  compute  a  posterior  probability  of  the  pixel  (n,m)  from  a  prior  probability  of  the 
same  pixel  and  a  set  of  prior  probabilities  of  the  adjacent  pixels.  The  prior  probabilities 
are  estimated  from  a  2-D  image  data  using  the  linear  filtering  model  (see  Equation 
2.25). 

The  Relaxation  formula  is  defmed  (Ref.  9]  by 


V  4- 1  Pr. 

Pjj  ^  ^  ^  (n,m)  =  Avg  (  — ■» - ^ - 


(A.I) 


where  k  is  the  number  of  the  iteration,  is  the  relaxation  factor,  P..*^  (n,m)  are  a  set 
of  prior  probabilities,  T  is  the  number  of  textures,  and  S  is  the  number  of  pixels.  The 
updated  estimates  P.j^  ^  (n,m)  are  obtained  by  averaging  all  of  the  terms  in 
parentheses.  The  relaxation  factor,  X,..*  ,  in  Equation  A.I  is  given  by 


(A.2) 


where  c(i,j|s,t)  is  a  nonnegative  compatibility  function  whose  value  is  small  if  the 
neighboring  pixel  is  black  when  the  estimated  pixel  is  white,  otherwise  its  value  will  be 
large.  X..*  is  used  to  update  the  probability  Pjj*^  (n,m).  Note  that  Xjj*  is  large  if  the 
compatibilities  c(i,j|s,t)  are  large  and  the  probabilities  are  high,  otherwise  Xy®  will 
be  small. 

The  results  of  the  Relaxation  Method  segmentation  are  shown  in  Figures  A.I 
and  A.2  .  The  Figure  A.I  presents  the  segmentation  of  Figure  4.7  with  the 
compatibility  c(i,x|s,x)  is  equal  to  0.1,  where  x  can  be  0  or  1.  The  result  is  very  spotty 
but  both  regions  are  discemable.  The  Figure  A.2  gives  the  segmentation  with  the 
compatibility  c(i,x|s,x)  is  equal  to  0.9.  This  result  is  better  than  the  previous  result,  but 
there  are  still  several  spots  especially  in  the  field  region.  Both  results  are  obtained  after 


41 


10  iterations.  I'igure  4.9  shows  the  result  of  MAP  segmentation  of  Figure  4.7  .  Ihe 
MAP  segmentation  result  is  much  more  accurate  than  the  Relaxation  method 
segmentation. 


Figure  A.l  Relaxation  Method  Segmentation  with  c  O.I. 


APPENDIX  B 

FILTER  PARAMETEI^|g|[J!^^<j[JJ^^D  MULTICHANNEL 


The  interactive  program,  FLTRl,  estimates  two  sets  of  filter  parameters  from  a 
2-D  three-channel  image.  The  corresponding  algorithms  were  presented  in  Section 
II.B.l.  In  order  to  run  this  program,  the  user  must  input  the  required  program 
parameters  in  the  order  listed  below 


le  number  of  filter  rows,  P.  Maximum  is  4. 

-he  number  of  filter  columns,  Q..  Maxunum  is  4. 

The  number  of  rows,  N,  m  each  unage  channel.  Maximum  is  128. 

The  number  of  columns,  M,  jn  each  unage  channel.  Maximum  is  128. 
The  number  of  channels,  K,  m  unage. 

The  filenarnes  of  the  image  chaimeis. 
he  coordinates  of  the  eslimation  windows  of  textures, 
wo  output  filenames  for  the  estimated  filter  parameters. 


The  mean  vectors  of  the  data  in  the  two  estimation  windows  is  first  computed  by 
FLTRl.  Then,  the  program  subtracts  the  mean  from  the  image  and  determines  the 
correlation  functions.  After  calculating  the  correlation  matrix,  the  program  determines 
the  covariance  matrix,  £  ,  using  the  equation 


I  R  ]  [  B  ]  »  [  S|  ]  (B  I) 

where  B  is  a  dummy  matra  that  has  the  same  dimension  with  A  matrix,  B^q  =  1 2^^  ]** 
and  Sj  is  all  zero  except  for  an  identity  matrix  in  its  first  partition.  The  covariance 
matrix  must  satisfy  the  relation 


B, 


00 


I 


Finally,  the  filter  weighting  coefficients  are  estimated  using 


(B.2) 


[A]-  [B].[Z] 


(B.3) 
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MULTICHANNEL  IMAGE  SEGMENTATION  PROGRAM 


The  interactive  program,  SGMTl,  segments  a  color  image  with  two  textures 
when  given  the  three-channel  image,  the  image  dimension,  and  two  sets  of  filter 


parameters.  Again,  the  user  must  input  the  required  program  parameters  to  run  the 


program.  These  parameters  have  to  be  in  the  order  listed  below  : 


le  number  of  rows,  N,  in  the  image  channels, 
le  number  of  columns,  M .  in  the  image  channels, 
le  number  of  filter  rows,  P. 
le  number  of  filter  colunms,  Q. 
le  number  of  textures,  T.  in  the  image, 
le  number  of  channels,  K,  in  the  image, 
le  filenames  of  the  ipiage  channels, 
le  filenames  of  the  filter  parameters, 
le  output  filename  of  the  ML  segmentation, 
le  output  filename  of  the  MAP  segmentation. 


This  program  estimates  the  errors  of  two  textures  using  Equation  2.7,  then  performs 


the  ML  segmentation  and  the  MAP  segmentation  using  Equations  2.25  and  2.33  .  The 
convergence  factor  KS,  and  the  size  of  must  be  assigned  properly  by  user  to 

perform  the  MAP  segmentation  accurately. 


APPENDIX  C 

KARHUNEN-LOEVE  TRANSFORMATION 


The  interactive  program,  KRLV,  implements  the  Karhunen-Lofeve  transformation 
algorithms  described  in  Section  III.B.  This  program  requires  the  user  to  assign  the 
program  parameters  in  the  order  listed  below  : 


each  channel. 

*  The  butput'fQenames  oFffic  tfansfonned  color  channels. 

The  correlation  matrix  is  first  calculated.  Then  the  program  calls  the  IMSL 
subroutine  EIGRS  to  calculate  the  eigenvalues  and  the  eigenvectors  of  the  correlation 
matrix.  Finally,  the  transformed  color  image  is  obtained  by  multiplying  the  transpose 
of  the  eigenvectors  matrix  by  the  original  color  image.  The  program  scales  the 
transformed  color  image  for  display  on  the  COMTAL  image  prossessing  system. 


ONE-CHANNEL  IMAGE  SEGMENTATION  PROGRAM 

The  program,  FLTR2,  calculates  two  sets  of  filter  parameters  from  a  2-D  single 
channel  image.  The  user  must  input  the  required  program  parameters  to  run  the 
program.  These  parameters  except  the  number  of  channel,  K,  are  given  in  Appendix  B. 

The  program,  SGMT2,  implements  a  single  channel  image  segmentation.  The 
required  program  parameters  given  in  Appendix  B  have  to  be  assigned  to  run  the 
program.  Equations  3.11  and  3.12  are  used  to  perform  the  ML  region  estimation  and 
the  MAP  region  estimation. 
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COMPUTER 


APPENDIX  D 


IMAGE 


A 

*  PROGRAM  FLTRl 

if 

A  PURPOSE  To  develop  two  sets  of  filter  parameters  from  a  2-D 

*  three-channel  input  image.  These  parameters  are  the 

*  mean  vector,  the  covariance,  and  the  filter  coeffici- 

*  ents  of  the  image. 

*  REQUIRED  IHSL  ROUTINES 

A  LEQT2F,  LUOATN,  LUELMN,  LUREFN 

*  IMPLEMENTED  BY  LTJG  TIMUR  KUPELI  Nov  1986 


aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

AAAA  VARIABLE  DEFINITIONS  **** 

BINPUT  3  [  Input  1  Image  data  in  byte  format. 

.FINPUT  s  f  Input  1  image  data  in  real*4  format. 

P  ,  Q  >  Rows,  Columns  of  the  linear  filter. 

N  ,  H  =  Rows,  Columns  of  the  image  data. 

K  >  The  number  of  channels  in  image. 

T  -  The  number  of  filters  for  processing. 

IMAGE  »  [  Input  ]  Filename  of  the  image  channels. 

FILTER  =  r  Output  ]  Filename  of  the  filter  parameters. 

MEAN  »  Array  of  mean  vectors  of  estimation  windows. 

R  The  correlation  matrix. 

IMATRX  =  Identity  matrix. 

KH  »  The  covariance  matrix 

A  =  The  filter  coefficients  matrix. 

ISIZE  >  The  estimation  window  size. 

NO,  NO  »  Row,  column  delay(shift)  of  correlation. 

This  program  uses  2  by  2  pixels  filter.  If  user  wants  to  use 
different  size  filter,  the  dimensions  of  the  following 
variables  must  be  modified. 

R  SI  A  WARE  A 

For  exempie.  ior'4  by  4  pixels  filter,  the  dimensions  must  be 
K(4d,48),  Slt48,3),  A(48,3},  WKAREA(2448) 


INTEGER*4  P , Q , ROW , COL , STARTN ( 2 ) , STARTM (2 ) , ENDN ( 2 ) , ENDM ( 2 ) , 

*  RNR0W,RNC0L,RR0W,CC0L,X,Y,K,1ER,QK,PKQ,L,J,JJ,JJJ, 

*  HNNO , HMMO , HCOL , HROW , LNNO , LNMO , LCOL , LROW , RN , RM , T , 

*  NO , MO , RRN , RRM , ROWl , COLl , MMO , NNO , IDGT , I 

REALM  SUN, TEMP, FINPUT (128, 128, 3) ,MEAN(2,3) ,WKAREA(180) , 

*  R(12,12),SI(12,3),IMATRX(3,3),B00(3,3),SUM1,SUM2,SUM3, 

*  KW(3,3),A(l2,35,fsiEE 

CHARACTER*50  IMAGE (3 ), FILTER (2) 

BYTE  BINPUT(128) 

T  »  2 


GET  PROGRAM  INPUT  PARAMETERS. 
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noon  o  o  non  o  non 


TYPE  10 

FORMAT ('  ENTER  NUMBER  OF  FILTER  ROWS  FROM  2-4  DESIRED  s' ,$) 
READ  11, P 
FORMAT (13) 

TYPE  12 

FORMATC  ENTER  NUMBER  OF  FILTER  COLUMNS  FROM  2-4  DESIRED 
READ  11, Q 

TYPE  13 

FORMATC  ENTER  NUMBER  OF  ROWS  IN  IMAGE  ' ,$) 

READ  11, N 

TYPE  14 

FORMATC  ENTER  NUMBER  OF  COLUMNS  IN  IMAGE  ',$} 

READ  11, M 

TYPE  15 

FORMAT?'  ENTER  NUMBER  OF  CHANNELS  [  MAX  -  3  ]  IN  IMAGE  ;',$) 
READ(*,11)  K 

GET  THE  MULTICHANNEL  IMAGE 
DO  16  J  =  1  ,  K 
WRITE(*,17)  J 

FORMATC  ENTER  FILENAME  OF  IMAGE  CHANNEL  '  13, $) 

READ(*,18)  IMAGE(J) 

FORMAT  (A50) 

CONVERT  THE  IMAGE  FROM  BYTE  FORMAT  TO  REAL*4  FORMAT 

OPEN(UNIT=l,FILE=IMAGE{J),STATUS='OLD' , ACCESS  =  'DIRECT') 

DO  180  ROW  s  1  ,  N 

READ  (I'ROW)  (BINPUT(COL),COL=l,M) 

DO  181  COL  =  1  ,  M 


CONTINUE 

CONTINUE 
CLOSE  (UNIT=1) 

CONTINUE 


TEMP  =  BINPUT(COL) 

IF  (TEMP.LT.O.Oj  THEN 
TEMP  =  TEMP  +  256 
END  IF 

FINPUT{ROW,COL,J)  =  TEMP 


GET  'T'  AREAS  FOR  WHICH  FILTERS  ARE  DESIRED  AND  OUTPUT  FILENAMES 
FOR  EACH  AREA'S  FILTER  COEFFICIENTS  ANi)  COVARIANCE  MATRIX. 

DO  19  I  =  1,T 
WRITE(*,2^  I 

FORMAT ( '  EIOTR  UPPER-LEFT  ROW  FOR  AREA  ' ,12, ' s ' ,$) 

READ<*,11)  STARTN(I) 

WRITE(*,21)  I 

FORMAT ( '  ENTER  UPPER-LEFT  COLUMN  FOR  AREA' ,12, ' s ' ,$) 
READ(*,11)  STARTM(I) 

WRITE(*,2^  I 

FORMAT ( '  ENTER  LOWER-RIGHT  ROW  FOR  AREA' ,12, ' i ' ,$) 

READ(*,11)  ENDN(I) 

WRITE (*,23)  I 

FORMATC  ENTER  LOWER-RIGTH  COLUMN  FOR  AREA' , 12, ' s ' , $) 
READ(*,11)  ENDM(I) 
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24 


26 

25 


WRITE (*,24)  I 


.$) 


FIND  THE  MEAN  VECTOR  OF  ESTIMATION  WINDOW  T  AREAS  OF  IMAGE 
CHANNELS.  THE  MEAN  VECTOR  CONSISTS  OF  THE  MEAN  FOUND  EACH  CHANNEL 


'^(ENDN(I)  -  STARTN(I)  +  1)*(ENDM(I)  -  STARTM(I)  +  1) 

o!o 
0.0 


ISIZE 
SUHl  3 
SUM2  « 

SUMS  » 

DO  25  L  »  STARTN(I),ENDN(I) 

DO  26  J  s  STARTM(I},ENOM(I) 
SUMl  *  SUMl  +  FINPUT(L,J,1 
SUM2  ■  SUM2  +  FINPUT(L,J,2 
*  FINPUTa.J,3 

CONTINUE 
CONTINUE 


HEAN(I,1 

MEAN(I,2 

MEAN(I.3 


SUMl  /  ISIZE 
SUM2  /  ISIZE 
SUM3  /  ISIZE 


—  suno  /  19J.6CI 

WRITE?*, 27)  (MEAN(I,II),II=1,K) 

27  FORMAT?'  MEAN* ' ,3F9.3,5j 

CORRECT  THE  IMAGE  TO  BE  ZERO  MEAN 

DO  2-^1  J  =  1  ,  K 

DO  28  L  =  STARTNd),  ENDN(I) 

DO  29  LL  f  STARTM(I),ENDM(I) 

29  CONTINUe’^'^^^'^^'^^  “  FINPUT(L,LL,J)  -  MEAN(I,J) 

28  CONTINUE 
271  CONTINUE 


291 


C 

C 


C 


DETERMINE  THE  2-D  CORRELATION  FUNCTION  OF  THE  IMAGE  CHANNEL 
TO  CORREIATION  MATRIX  APPROPRIATE  TO  THE  INPUT  INPUTpStERS 

UP  r,  R,  g 


CORRELATION . MATRIX ’ , I 2 , $ ) 


WRITE(*,291)  I 
FORMAT ( ’ 

PKQ  s  P  *  K  *  Q 
OK  =  Q  *  K 
DO  30  RNROW  »  1,P 
X  »  OK  *  (RNROW-1) 

RNCOL  =  1,P 
NO  =  RNROW  -  RNCOL 
Y  *  QK  *  <RNCOL-l) 

DO  32  RMROW  *  1,0 

RN  »  K  *  (RMROW  -  1)  +  X 
DO  33  RMCOL  =  1,0 

MO  =  RMROW  -  RMCOL 
RM  =  K  *  (RMCOL  -  1 


DO  3 


)  +  Y 


LNNO  =  STARTN(I)  +  NO 
LMMO  =  STARTM(I)  + 


HNNO  =  ENDN 
HMMO  =  ENDM 


LCOL  -I 


ai: 


NO 

MO 


MO 


HCOL 

wow  =MAXO(STARfN(I),LNNO) 
HROW  =MINO(EMDN(I),HNNO) 


DO  133  R0W1=  1  ,  3 
RRN  »  RN  +  ROWl 

DO  233  COLl  =1,3 
RRM  =  RM  +  COLl 
SUM  «  0.0 

DO  333  ROW  =  LROW  ,  HROW 
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NNO  =  ROW  -  NO 

IF  (NNO  .LT.  LROW)  THEN 

RROH  s  HAXO(ROH,NNO) 
ELSE  IF  (NNO  .GT.  HROW)  THEN 
RROH  *  MIN (ROW, NNO) 


ELSE 
END  IF 
DO  433 


C 

C 


RROH  =  NNO 


COL  =  LCOL  ,  HCOL 
MMO  =  COL  -  MO 
IF  (MMO  .LT.  LCOL)  THEN 
CCOL  »  MAXO(MMO,COL) 

ELSE  IF  (MMO  .GT.  HCOL)  THEN 
CCOL  =  MINO(MMO,COL) 

ELSE 

CCOL  =  MMO 

END  IF 


SUM  =  SUM  +  FINPUT(ROW,COL,ROWl)  *  FINPUTSAVEW,CCOL,COLl) 


433 

333 

233 

133 

33 

32 

31 

30 


SUM  /  ISIZE 


CONTINUE 
CONTINUE 
R(RRN,RRM) 

CONTINUE 
CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

THE  FOLLOWING  FORMAT  MUST  BE  MODIFIED  TO  USE  DIFFERENT  FILTER  SIZE 
THEN  2  by  2  PIXELS. 

DO  330  L  =  1  ,  PKQ 
WRITE?*, 331)  (R(L,J),J=1,,PKQ) 


331 

FORMAT?' 

332 

FORMAT?' 

330 

CONTINUE 

12F6.0) 


37 

36 


RESET  THE  IDENTITY  MATRIX. 

DO  36  J  =  1  ,  K 
DO  37  L  =  1  ,  K 
IF  (J.EQ^.L) 

IMATRX{ 

ELSE 

IMATRX(J,L)  =  0.0 
END  IF 

CONTINUE 

CONTINUE 


L)  THEN 
C(J,L)  =  1.0 


RESET  SI(J,L)  TO  HAVE  [  I  ]  IN  FIRST  PARTITION  AND  [  0  ]  IN  ALL  OTHERS 

DO  38  J  =  1  ,  PKQ 
DO  39  L  =  1  ,  K 

IF  (J.EQ.L)  THEN 
SI(/,L)  =1.0 
ELSE 

SI(J,L)  =  0.0 
END  IF 

39  CONTINUE 

38  CONTINUE 


SOJ-y^-IQy^TION^^R^^  lI'.Ll-l 


CALL  TO  IMSL  ROC 


I  SI  ]  .  NOTE  THAT  THE  BELOW 


LEQT2F,  (  B  ]  RETURNS  IN  [  SI  ] 
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CALL  LEQT2F  (R,K,PKQ,PKQ,SI,IDG,WKAREA,IER) 

SOLVE  EQUATION  [  BOO  1  *  f  KW  1  *  [  I  1  FOR  COVARIANCE 
MATRIX  [  KW  1  AFTER  COLLING  IMSL  ROUTINE  LEQT2F,  [  KW  1 
RETURNS  IN  [  IMATRX  ].  y  /  l  J 

DO  45  J  =  1  ,  K 

DO  46  JJ  s  1  ,  K 

A/;  “  SI{J,JJ) 

46  CONTINUE 
45  CONTINUE 

CALL  LEQT2F ( BOO , K , K , K , IMATRX , IDG , WKAREA , lER ) 

SOLVE  FOR  THE  FILTER  COEFFICIENTS  ,  [A]=[B]*[KWl, 

WHICH  IN  THE  PROGRAM  IS  [  A  ]  =  [  SI  ]  *  [  iNaTRX  ] 

DO  47  J  =  1  ,  PKQ 
DO  48  JJ  *  1  ,  K 
TEMP  =0.0 
DO  49  JJJ  =  1  ,  K- 

TEMP  =  TEMP  +  SI(J,JJJ)  *  IMATRX(JJJ,JJ) 

49  CONTINUE  ^  ' 

„  A(J,JJ)  =  TEMP 

48  CONTINUE 

47  CONTINUE 

WRITE/*, 491)  ((IMATRX(J,JJ),JJ=1,K),J=1,K) 

491  F^TC  •,  <K>(F6.2,3X)) 

WRITE{*,53)  I 

53  FORMATC  FILTER  COEFFICIENTS  12, $) 

WRITE(*,491)  ((A(J,JJ),JJ=1,K),J=1,PKQ) 

WRITE  OUT  MEAN,  COVARIANCE  MATRIX,  AND  FILTER  COEFFICIENTS  TO  THE 
USER  INPUT  FILE. 

OPEN  OTJ=2 ,STATUS='NEW'  , CARRIAGECONTROL= ' LIST '  , 
*  FORM= ' FORMATTED ' ) 

Aoc  (MEAN(I,J),J=1,K)  • 

495  FORMAT (FIO. 4) 

WRITE(2,495)  ( (IMATRX(J, JJ) , JJ=1 ,K) , J=1 ,K) 

WRITE(2,495)  aA(J, JJ) , JJ=i ;k) , J=1,PKQ) 

’  CLOSE (UNIT=2) 

19  CONTINUE 

WRITE(*,55) 

55  FORMATC  THE  PROGRAM  FLTRl  IS  OVER',$) 

STOP 

END 
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****  VARIABLE  DEFINITIONS  **** 

BINPUT  =  Input  ]  Image  data  in  the  byte  format. 

ML  =  Output  ]  The  result  of  ML  segmentation  in  byte  format. 

MAP  =  Output  J  The  result  of  MAP  segmentation  in  byte  format. 

FNAME  =  '  Input  ]  Filename  of  filter  parameters  set. 

IMAGE  =  Filename  of  the  image  channel. 

P  ,  Q  =  Rows  ,  Columns  of  filter. 

N  ,  M  =  Rows  ,  Columns  of  image. 

K  =  Number  of  channel  of  image. 

KW  =  [  Input  1  The  covariance  matrix. 

MEAN  =  [  Input  J  Mean  vector  of  estimation  windows. 

ERROR  =  Prediction  error  estimation 
A  =  [  Input  }  Filter  coefficients  matrix. 

TEXTURE  =  Zero-mean  image  data*  in  real*4  format. 

COl  =  Number  of  removed  false  points  in  the  first  texture. 

CIO  =  Number  of  removed  false  points  in  the  other  texture. 

IN  ,  IM  =  Maximum  number  of  rows  and  columns  in  the  image. 

PI  ,  Q1  =  Maximum  number  of  rows  and  columns  in  the  filter. 

TMAX  =  Maximum  number  of  filters  for  processing. 

IK  =  Maximum  number  of  channels  in  image. 


1 

2 

3 

C 


INTEGER  IN,N,IM,P1,P^1,Q,TMAX,T,PQ,R0W,C0L,I,J,JJ,JJJ,L,LL,KK, 
LLL , LLLL . COUNT; LI ,HI ,LJ ,HJ , K, PP , QQ , IK , II , III , COl , CIO , M 

REAL  KW(l:3,ls3,l;2),MEAN<l.3,l»2),TEMP,SUMl,SUM2,PMLll,LN(2), 
ERR0R(1 : 128,  ia28,ls3,l:2),PMLl,PNL2,PML(l:  128, 1:128), 
AREA,KS,A(ls3,ls3,l-.2,l:2,l-.2),Dl,D2,KWl(l*.3,l:3),KW2(3,3), 
EKWl(l,l:3),EKW27l,l:3),PML22,DD2,TEXTURUtl28,l:128,3,3) 
WKA^A(6),Aa(1:128;1:128) 

CHARACTER*50  IMAGE (1:3), FNAME (1:2), MLTEST , MAPTEST 

BYTE  BINPUTU :128, 1 :128, 1;3),ML(1;128,1 -.128)  ,MAP(1:128, 1:128) 
MLI(1 :l28, 1 :128) 

IN  =  128 
IM  =  128 
PI  =  4 
01  =  4 
TMAX  3  2 
IK  =  3 

GET  THE  INPUT  PARAMETERS  OF  THE  PROGRAM 


WRITE(*,2)  IN 

FORMAT ('  ENTER  THE  NUMBER  OF  ROWS  IN  IMAGE. LIMIT  OF ' , 13 , '  : ' , $) 
READ(*,3)  N 
F0RMAT(I3) 

IF((N.LT.l)  .OR.  (N.GT.IN))  GOTO  1 
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WRITE(*,5)  IM 

FORMATC  ENTER  THE  NUMBER  OF  COLUMNS  IN  IMAGE. LIMIT  0FM3, 
READ(*,3)  M 

IF((M.LT.l)  .OR.  (M.GT.IM))  GOTO  4 
WRITE(*,7)  PI 

FORMATC*  ENTER  THE  NUMBER  OF  ROWS  IN  FILTER. LIMIT  0FM3, 
READ(*,3)  P 

IF((P.LT.2)  .OR.  (P.GT.Pl))  GOTO  6 
WRITE(*,9)  01 

FORMATC  ENTER  THE  NUMBER  OF  COLUMNS  IN  FILTER. LIMIT  0FM3, 
READ(*,3)  Q 

IF((Q.LT.2)  .OR.  (Q.GT.Ql))  GOTO  8 
WRITE(*,11)  TMAX 

FORMATC  ENTER  NUMBER  OF  TEXTURES  TO  PROCESS .LIMIT  OF ' , 13 , ' : ' , 
READ(*,3)  T 

IF((T.LT.2)  .OR.  (T.GT.TMAX))  GO  TO  10 
WRITE(*,13)  IK 

FORMATC  ENTER  THE  NUMBER  OF  IMAGE  CHANNELS. LIMIT  OF','  ',I3,$) 
READ(*,3)  K 

IF((K.LT.l)  .OR.  (K.GT.IK))  GO  TO  12 

GET  ALL  CHANNELS  OF  THE  IMAGE 

DO  19  I  =  1  ,  K 
WRITE(*,20)  I 

FORMATC  ENTER  FILENAME  OF  THE  IMAGE  CHANNEL  NUMBER','  ',I2,$) 
READ(*,21)  IMAGE(I) 

FORMAT (A50) 

OPEN(UNIT^l,FILE=IMAGE(I) ,STATUS='OLD' ,ACCESS= ' DIRECT ' ) 

DO  23  ROW  =  1  .  N 

READ(l'ROW)  (BINPUT(ROW,COL,I),COL  =  1  ,  M  ) 

CONTINUE 

CLOSE (  UNIT  =  1  ) 

CONTINUE 

GET  THE  FILTER  COEFFICIENT, MEAN,  AND  COVARIANCE  MATRICES 

DO  25  1=1, T 
WRITE(*,26)  I 

FORMATC'ENTER  FILENAME  OF  FILTER  PARAMETERS  SET  NUMBER','  ',I2,$) 
READ(*,21)  FNAME(I) 

OPEN(UNIT=2 , FILE=FNAME ( I ) , STATUS® ' OLD ' , FORM= ' FORMATTED ' ) 

DO  260  J  =  1  ,  K 
READ(2,27)  MEAN(J,I) 

FORMAT (FIO. 4) 

CONTINUE 

DO  28  J=1,K 

READ(2,27)  (KW(J, JJ,I) , JJ=1,K) 

CONTINUE 

DO  29  PP=1,P 
DO  30  QQ=1,Q 
DO  31  ROW=l,K 

READ (2,33)  ( A ( ROW , COL , PP , QQ , I ) , COL=l , K ) 

FORMAT (FIO. 4) 

CONTINUE 

CONTINUE 

CONTINUE 
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32 

C 

25 


220 


37 

36 

35 


46 

45 

44 

43 

421 

42 

41 

40 

C 


48 

47 

C 


C 


CLOSE (UNIT=2) 

FORMAT{3(F10.4,3X)) 

CONTINUE 

GET  THE  OUTPUT  FILENAMES  OF  THE  ML  AND  HAP  SEGMENTATION  RESULTS. 

READ (*,220)  MLTEST 
READ (*,220)  MAPTEST 
FORMAT (A80) 

CONVERT  A  2-D  BYTE  INPUT  IMAGE  DATA  ARRAY  IN  THE  RANGE  OF  -128  TO 
127  THAT  REPRESENTS  APPROPRIATE  INTENSITY  LEVELS  IN  THE  RANGE  OF 
0  TO  255  . 

DO  35  J  =  1  ,  K 
DO  36  ROW  =  1  ,  N 

DO  37  COL  =  1  ,  M 

TEMP  =  BINPUT(ROW,COL,J) 

IF(TEMP.LT.O.O)  THEN 

TEMP  =  TEMP+256 

END  IF 

TEXTUR(ROW,COL,J,l)  =  TEMP  -  MEAN(J,1) 
TEXTUR(R0W,C0L,J,2)  =  TEMP  -  MEAN(J,2) 

CONTINUE 

.  CONTINUE 
CONTINUE 

CALCULATION  OF  ERROR  ESTIMATE  FOR  TWO  TEXTURES 

DO  40  I  =  1  ,  T 

DO  41  L  =  1  ,  N 
DO  42  LL  =  1  ,  M 
♦  DO  421  KK  =  1  ,  K 

ERROR(L,LL,KK,I)  =  0.0 
DO  43  III  =  1  ,  P 
J  *  1  -  III 
DO  44  LLL  =  1  ,  Q 
JJJ  a  LL  -  LLL 
DO  45  II  a  1  ,  K 
DO  46  JJ  =  1  ,  K 
IF(J.LE.O)  Jal 
IF(JJJ.LE.O)  JJJ=1 

ERROR(L,LL,KK,I)=ERROR(L,LL,KK,I)+A(II,JJ,III,LLL,I)*TEXTUR(J,JJJ,KK,I) 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

DO  47  JJ=1,K 
DO  48  LL=1,K 

KW1(JJ,LL)  =  KW(JJ,LL,1) 

KW2(JJ,LL)  =  KW(JJ,LL,2) 

CONTINUE 

CONTINUE 

D1  =  1.0 

CALL  LINV3F(KW1 ,6, 1 ,K,K,D1 ,D2,WKAREA, lER) 

DETl  a  D1  *  2**D2 
LNU)  =  ALOG(DETl) 

D1  a  1.0 

CALL  LINV3F(KW2,6,1,K,K,D1,DD2,WKAREA,IER) 
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DET2  =  D1  *  2**DD2 
LN(2)  =  AL0G(DET2) 

CALCULATION  OF  MAXIMUM-LIKELIHOOD  IMAGE  SEGMENTATION 

OPEN (UNIT=3 , FILE=MLTEST , STATUS* ' NEW ' , ACCESS* ' DIRECT ' , 

*  RECL*(IM/4),MAXREC*IN) 

DO  50  ROW  *1,N 
DO  51  COL  *  l.M 
DO  52  I  *  1,K 

EKW1(1,I)  =  0.0 
EKW2a,n  =  0.0 
DO  53  J*1,K 

EKWl  (1,1) *EKW1 (1,1) +ERR0R (ROW , COL , J , 1 ) *KH1 
EKW2 a , I } *EKH2 a , I ) +ERROR ( ROW , COL , J , 2 ) *RW2 

53  CONTINUE 

52  CONTINUE 

PMLll  =  0.0 
PML22  *  0.0 
DO  54  L*1,K 

PMLl 1=PML1 1+EKWl ( 1 , L ) *ERROR ( ROW , COL , L , 1 ) 
PML22*PML22+EKW2  U . L) *ERR0R ( ROW , COL , L , 2 ) 

54  CONTINUE 

PMLl  =0.0 
PML2  =0.0 
PMLl  =  PMLll  +  LN(1) 

PML2  =  PML22  +  LNU) 

ML(ROW,COL)=0 
PML(ROW,COL)  =  PML2  -  PMLl 

IF(PML1  .GT.  PML2)  THEN 
fa(ROW,COL)=  -1 
END  IF 

MLI(ROW,COL)  *  ML(ROW,COL) 

51  CONTINUE 

WRITEO'ROW)  (ML(ROW,COL),COL*l,M) 

50  CONTINUE 

CLOSE (UNIT=3) 

MAXIMUM  A  POSTERIORI  IMAGE  SEGMENTATION 

0PEN(UNIT=4 , FILE«MAPTEST , STATUS= ' NEW ' , ACCESS* ' DIRECT ' 

*  RECL=  (IM/4),  MAXREC  *  IN) 

KS  »  10.0 
COl  *  0 
CIO  *  0 

DO  600  II  *  1  ,  5 
DO  60  1=1, N 

DO  61  J  =  1  ,  M 
SUMl  =  0.0 
SUH2  =0.0 
LI  =  I  -  3 
HI  *  I  +  3 
LJ  =  J  -  3 
HJ  *  J  +  3 

IF(LI.EQ.O)  LI  =  1 
IF(HI.GT.N)  HI  =  N 
IF(LJ.EQ.O)  LJ  =  1 
IF(HJ.Gt.M)  HJ  =  M 

AREA  =  (HI  -  LI  +  1)  *  (HJ  -  LJ  +  1) 

DO  62  ROW  =  LI  ,  HI 
DO  63  COL  =  LJ  ,  HJ 
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62 


CONTINUE  ■  MI-KROW.COL) 

CONTINUE 

END 


61  CONTINUE 
60  CONTINUE 
600  CONTINUE 


END  IF 
Ml-Id.J)  =  MAP(I,J) 


71 

70 

64 


DO  JO  I  »  1  ,  N 

DO  71  J  =  'i  ,  M 

ELSE  •^*  «D(I,J))  COl  = 

^IF<JIM(I,J)  .HE.  m.(l 

,J))  C10=C10+1 

CONTINUE 

CONTIITOE  (l<»P(I,J).J-l.M) 

CLOSE (UNIT=4) 

coi®?yi^*§^t)ciO:{?is) 

END 


COl  +  1 
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APPENDIX  E 

PROGRAMS  FOR  K-L  TRANSFORMATION,  AND  ONE-CHANNEL  SEGMENTATION 
c 

Q  ic  * 

C  *  PROGRAM  KRLV  * 


*  PURPOSJE  To  impl.ement  Karhunen-Loeve  transformation  from  a  2-D 

*  color  Image  which  size  is  128  by  128  pixels.  * 

*  REQUIRED  IMSL  ROUTINES  * 

*  EIGRS,  EQRT2S,  EHOBKS,  EHOUSS,  UERTST,  USPKD,  UGETIO  * 

*  IMPLEMENTED  BY  LTJG  TIMUR  KUPELI  Dec  1986  * 

A  A 

INTEGER  I,J,L,K,N1,N2,NN,R0W,C0L,ITEMP,QC(128,128), 

*  MAXVAL,MINVAL,IDIF,INTVAL,N 

REAL  R(l:3,l:3),E(ls3,l:3),FINPUT(l:128, 1:128, 1:3), 

*  D(3), WK(3), TEMP, Q{1;128, 1:128, 1:3), SLOPE 

CHARACTER*50  IMAGE(1:3)  ,  FNAME(1:3) 

BYTE  BINPUT(1:128,1:128),QQ(1:128, 1:128, 1:3) 

TYPE  100 

FORMAT (  'ENTER  THE  NUMBER  OF  ROWS,  COLUMNS  IN  THE  IMAGE ',$) 

READ  101 ,M 
FORMAT (13) 

TYPE  102 

FORMAT ('  ENTER  THE  NUMBER  OF  CHANNELS  IN  THE  IMAGE ',$) 

READ  101, K 

GET  THE  FILENAMES  OF  THE  RED, GREEN, AND  BLUE  COMPONENTS 
OF  THE  COLOR  IMAGE. 


DO  1  I  =  1  ,  K 
WRITE?*, 2)  I 

FORMAT ('  ENTER  THE  FILENAME  OF  THE  IMAGE  CHANNEL  NUMBER  ' 


WRITE?*, 2)  I 
ORMAT?'  ENTER  THE  FI 
READ(^,3)  IMAGE (I) 
WRITE?*, 3)  IMAGE(l) 
WRITE (*.45) 
ORMAT(A50) 


,I2,$) 


WRITE (*,45) 

FORMAT (A50) 

CONVERT  THE  IMAGE  BTYE  FORMAT  TO  THE  REAL  NUMBER  FORMAT 
OPEN(UNIT=l,  FILE=IMAGE(I),  STATUS='OLD' ,  ACCESS* ' DIRECT ' ) 
DO  5  ROW  *  1  ,  N 

READ  ( 1 '  ROW  )  ( BINPUT  ( ROW ,  COL  )  ,  COL=l ,  N) 

DO  6  COL  *  1  ,  N 

TEMP  =  BINPUT (ROW, COL) 

IF  (TEMP.LT.O.O)  THEN 

TEMP  =  TEMP  +  256 

END  IF 

FINPUT(ROW,COL,I)  =  TEMP 
CONTINUE 

CONTINUE 
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CLOSE (UNIT  =  1) 


CONTINUE 

CALCULATE  THE  CORRELATION  MATRIX 

NN  =  N  *  N 
DO  20  I  =  1  ,  K 
DO  21  J  =  1  ,  K 
R{I,J)  =  0 
DO  22  N1  =  1  ,  N 
DO  23  N2  =  1  ,  N 

R(I,J)  =  R(I,J)  +  FINPUT(N1,N2,I)*FINPUT(N1,N2,J) 
CONTINUE 
CONTINUE 

R(I,J)  =  R(I,J)  /  NN 
CONTINUE 
CONTINUE 

WRITE(*,45) 

WRITE (*.30) 

WRITE(*,39) 

FORMATC  . ') 

FORMAT ('  ' ,5X, '  THE  CORRELATION  MATRIX' ,$) 

WRITE/*, 31)  ((R(I.J),Jsl,K),I=l,K) 

FORMAT(  <K>(F9.2,4X)) 

WRITE (*,39) 

CALCULATE  EIGENVALUES  AND  EIGENVECTORS  OF  THE 
CORRELATION  MATRIX 

JOBN  =11 

CALL  EIGRS(R,K,JOBN,D,E.K,WK,IER) 

SORT  THE  EIGENVALUES  IN  DECREASING  ORDER 


TEMP  =  D(l) 
dB!  *  TEMP 


DO  49  I  =  1  .  K 

TEMP  =  E(I,1) 
e7i,1)  =  Eh, 3) 
E(i,3)  =  TEMP 

CONTINUE 

WRITE(*,39) 

WRITE (*,45) 

FORMATC  ') 


WRITE (*.41) 

WRITE (*, 39) 

FORMATC  ',5X,'  THE  EIGENVALUES  ',$) 
DO  46  J=1,K 
WRITE(*,42)  D(J) 

FORMAT (5X,F9. 2) 

CONTINUE 
WRITE (*,39) 

WRITE (*,45) 

WRITE (*,43) 

WRITE (*,39) 

FORMATC  ',5X,'  THE  EIGENVECTORS',  $) 
WRITE/ *.44)  ((E(I,J),J=1,K),I=1,K) 
FORMAT(3(F9,2,4X)) 

WRITE (*,39) 


THE  EIGENVALUES  ',$) 


USE  THE  TRANSPOSE  OF  THE  EIGENVECTORS  MATRIX  TO  IMPLEMENT 
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onn  non  o  no 


KARHUNEN-LOEVE  TRANSFORMATION. 


DO  50  I  *  1  ,  K 

DO  51  N1  =  1  ,  N 

DO  52  N2  =  1  ,  N 
g(Nl^N2,I)  =  0.0 

DO  53  J  =  1  ,  K 

Q(N1,N2,I)  =  Q(N1,N2,I)  +  E(J,I)*FINPUT(N1,N2,J) 
CONTINUE 
CONTINUE 
CONTINUE 
CONTINUE 

DO  60  I  a  1  ,  K 
WRITE(*,61)  I 

FORHATC  ENTER  FILENAME  OF  THE  TRANSFORMED  IMAGE  ',I2,$) 

READ(*,3)  FNAMEjfl) 

WRITER*, 3^  FNAME(I) 


WRITE 


CONVERT  THE  TRANSFORMED  IMAGE  TO  BYTE  FORMAT 

DO  62  N1  =  1  ,  N 
DO  63  N2  =  1  ,  N 

QC(N1,N2)  =  JNINT(Q(N1,N2,I)) 

CONTINUE 

CONTINUE 

SCALE  THE  TRANSFORMED  IMAGE  TO  BE  WITHIN  DISPLAY  RANGE 

IF  (I  .EQ.  1)  THEN 
MAXVAL  =  QC(1,1) 

MINVAL  =  QCU,1) 

DO  ROW  =  1  ,  N 
DO  COL  =  1  ,  N 

IF  (QC(ROW,COL)  .GT.  MAXVAL)  THEN 
fftXVAL  =  QC (ROW, COL) 

ELSE  IF  (QC<ROW.COL)  .LT.  MINVAL)  THEN 
MINVAL  a  QC (ROW, COL) 

END  IF 
ENDDO 
ENDDO 

INTVAL  a  MAXVAL  -  MINVAL 
SLOPE  a  255  /  REAL (INTVAL) 

DO  ROW  a  1  ,  N 
DO  COL  a  1  ,  N 

IF  (QC (ROW, COL)  .EQ.  MINVAL  )  THEN 
QC(ROW,COL)  a  0 

ELSE  IF  (QC(ROW,COL)  .EQ.  MAXVAL)  THEN 
QC(ROW,COL)  a  255 

pr  gg 

IDIF  a  QC(ROW,COL)  -  MINVAL 
QC(ROW,COL)  a  INT(SLOPE*IDIF) 

END  IF 
ENDDO 
ENDDO 
ELSE 

MINVAL  a  QC(1,1) 

DO  ROW  a  I  ,  N 
DO  COL  a  1  ,  N 

IF  (QC(ROW,COL)  .LT.  MINVAL)  THEN 
mNVAL  a  QC  (ROW,  COL) 

END  IF 
ENDDO 
ENDDO 

DO  ROW  a  1  ,  N 
DO  COL  a  1  ,  N 

IDIF  a  QC(ROW,COL)  -  MINVAL 


r 


oo 
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64 

C 


ENDOO 
ENDDO 
END  IF 


QC{R0W,C0L) 


DO  64  N1 
DO  65  N2 


CONTINUE 

CONTINUE 


=  1  ,  N 

END  IP  -  "Elf 

QQ(m,N2,I)  =  ITEM? 


66 

C 

C 

60 

C 


DO  66  N1  s  I  M 

(QQ(N1,N2,I),N2=1,N) 


CLOSE (UNIT®  2) 

CONTINUE 

STOP 

END 


60 


t 


1NT{SL0PB*IDIF) 


-  256 


NEW .ACCESS® 'DIRECT' , 


non  ono  o  n  o  onn 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 


C 

C 

C 


A 

*  PROGRAM  FLTR2 

A 

*  PURPOSE  To  develop  two  sets  of  filter  parameters  from  a  2-D  single 

*  channel  image  which  size  is  12B  by  128  pixels.  These  para- 

*  meters  are  the  mean,  the  covariance,  and  the  filter 

*  coefficients. 

•k 

*  REQUIRED  IMSL  ROUTINES 

A  LEQT2F,  LUDATN,  LUELMN,  LUREFN 

*  IMPLEMENTED  B7  LTJG  TIMUR  KUPELI  Sep  1986 
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 


A 

A 


A 

A 


INTEGER*4  P , Q , ROW , COL , STARTN ( 2 ) , STARTM ( 2 ) , ENDN ( 2 ) , ENDM ( 2 ) , IDGT , 

PQ , RNROW , RNCOL , RROW , CCOL , X , Y , 

HNNO ,HMm6 ,HC0L ,HROW,LNNO , LNMO , LCOL , LROW , RN , RM , lER 


REAL*4  SUM,TEMP,FINPUT(128,128) ,MEAN(2) ,WKAREA(28) , 
R(4.4),SI(4,1),IMATRX{1,1),B00, 
KW(i,lJ,AU'l),ISI2E 

CHARACTER*50  IMAGE ,FILTER(2) 

BYTE  BIMPUT(128) 


T  »  2 


GET  PROGRAM  INPUT  PARAMETERS*. 

TYPE  10 

10  FORMATC  ENTER  NUMBER  OF  FILTER  ROWS  FROM  2-4  DESIRED  «',$) 

READ  11. P 

11  F0RMAT(I3) 

TYPE  12 

12  FORMATC  enter  NUMBER  OF  FILTER  COLUMNS  FROM  2-4  DESIRED 

READ  11, Q 

TYPE  13 

13  FORMATC  ENTER  NUMBER  OF  ROWS  IN  IMAGE  ',$) 

READ  11, N 

TYPE  14 

14  FORMATC  ENTER  NUMBER  OF  COLUMNS  IN  IMAGE  ' ,$) 

READ  11, M 

GET  FILENAME  OF  SINGLE  CHANNEL  IMAGE 
TYPE  15 

15  FORMATC  ENTER  FILENAME  OF  IMAGE  ',$) 

READ  16, IMAGE 

16  FORMAT (A50) 

CONVERT  THE  IMAGE  BTYE  FORMAT  TO  REAL  NUMBER  FORMAT 

OPEN(UNIT=l , FILE»IMAGE , STATUS* ' OLD ' , ACCESS  *  ' DIRECT ' > 

DO  17  ROW  *  1  ,  N 

READ  (I'ROW)  {BINPUT(COL),COL=l,M)  ’ 

DO  18  COL  »  1  ,  M 

TEMP  =  BINPUT{COL) 

IF  (TEMP.LT.O.Oj  THEN 
TEMP  *  TEMP  +  256 
END  IF 

FINPUT(ROW,COL)  *  TEMP 
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18 

17 


20 


21 


22 


23 


24 


26 

25 


CONTINUE 

CONTINUE 
CLOSE  (UNIT=1) 

GET  'T'  AREAS  FOR  WHICH  FILTERS  ARE  DESIRED  AND  OUTPUT  FILENAME  FOR 
EACH  AREA'S  FILTER  COEFFICIENTS  AND  COVARIANCE  MATRIX. 

DO  19  I  =  1,T 
WRITE (*.20)  I 

FORMAT ('  ENTER  UPPER-LEFT  ROW  FOR  AREA  ',12,':'.$) 

READ(*,11)  STARTN(I) 

WRITE (*,21)  I 

FORMAT ('  ENTER  UPPER-LEFT  COLUMN  FOR  AREA' ,12, ' i ' ,$) 

READ(*,11)  STARTM(I) 

WRITE(*,22)  I 

FORMATC  ENTER  LOWER-RIGHT  ROW  FOR  AREA' ,12, ' i ' ,$> 

REAO(*.ll)  ENDN(I) 

WRITE(*,23)  I 

FORMAT ( '  ENTER  LOWER-RIGTH  COLUMN  FOR  AREA' ,12, ' : ' ,$) 

READ(*,11)  ENDM(I) 

WRITE(*,24)  I 

FORMAT ( '  ENTER  OUTPUT  FILE-SPEC  FOR  FILTER' ,12, ' ; ' ,$) 

READ (*,16)  FILTER(I) 

FIND  THE  MEAN  VECTOR  OF  ESTIMATION  WINDOW  AREA  OF  IMAGE 

ISIZE  =  ^ENDN(I)  -  STARTN(I)  +  1)*(ENDM(I)  -  STARTM(I)  + 

DO  25  ‘l  =  STARTN(I),ENDN(I) 

DO  26  J  a  STARTM(I),ENDM(I) 

SUM  a  SUM  +  FIOTUT(X,J) 

CONTINUE 
CONTINUE 

MEAN  (I)  a  SUM  /  ISIZE 
WRITE(*,27)  I.MEANd) 


1) 


U  t  I  *  «  \  ^  / 

27  FORMATC  MEANC '(' .12, ')',':' ,F9. 2, $) 

CORRECT  THE  IMAGE  TO  BE  ZERO  MEAN 


29 

28 


DO  28  La  STARTN(I),  ENDN(I) 

DO  29  J  a  STARTM(1),ENDM(I) 

FINPUT(L,J)  - 

CONTINUE 
CONTINUE 


FINPUT(L,J)  -  MEAN(I) 


291 


DETERMINE  CORRELATION  MATRIX 
WRITE(*,291)  I 

FORMAT ( '  CORRELATION . MATRIX ' , I 2 , $ ) 


a  P  * 
30  rn: 
X  a 
DO  3 


How  a  i,p 
■  (RNROW-1) 

RNCOL  a  1,P 
NO  a  rnROW  -  RNCOL 


Y  a 
DO  3 


*  (RNCOL-1) 

RMROW  «  1,0 
RN  a  X  +  RMROW 
DO  33  RMCOL  =1,0 

NO  a  RMROW  -  RMCOL 
RM  a  Y  +  RMCOL 
LNNO  a  STARTN 
LMMO  a  STARTM 
HNNO  a  eNDN 
HMMO  a  eNDM 


a! : 


lai: 


NO 

MO 


MO 
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LCOL  =1 
HCOL  =1 
LROW  =1 
HROW  =1 


=KAX0 ( STARTM ( I ) , LMMO ) 
=MiN0  ( ENDM  ( I ) ,  mmo ) 
=MAXO ( STARTN ( I) , LNNO ) 
=MIN0(ENDN(I),H^0) 


SUM  =  0.0 

DO  34  ROW  =  LROW, HROW 

NNO  =  ROW  -  NO 
DO  35  COL  =  LCOL, HCOL 
MMO  =  COL  -  MO 

SUM  =  SUM  +  FINPUT(ROW,COL)*FINPUT(NNO,MMO) 

35  CONTINUE 

34  CONTINUE 

R(RN,RM)  =  SUM  /  ISIZE 

33  CONTINUE 

32  CONTINUE 

31  CONTINUE 

30  CONTINUE 

RESET  SI(J,1)  TO  HAVE  1  IN  THE  FIRST  ROW  AND  0  IN  ALL  OTHERS 

DO  41  J  =  1,PQ 

IF  (J.EQ.l)  THEN 

SI(J,1)  =  1.0 
IMATRXU/1)  =  1*0 

ELSE 

SI(J,1)  =  0.0 

END  IF 

41  CONTINUE 

THE  FORMATS  BELOW  MUST  BE  MODIFIED  TO  USE  DIFFERENT  FILTER 
SIZE  THAN  2  by  2  PIXELS. 

IDGT  =1 
DO  36  K=1,PQ 


WRITE 

37 

FORMAT ( ' 

38 

FORMAT r 

36 

CONTINUE 

:•  4F9.2,! 


SOLVE  EQUATION  [  R  ]  *  [  B  ]  =  [  SI  ] 

CALL  LEQT2F  (R, IDGT,PQ,PQ,SI ,3,WKAREA, lER) 

WRITE(*,360)  (SI(J,1),J=1,PQ) 

FORMAT (FIO. 2) 

SOLVE  EQUATION  BOO  *  KW  =  I 
BOO  =  SI(1,1) 

GET  THE  FILTER  COEFFICIENTS  USING  THE  EQUATION 
t  A  ]=  [  B  ]  *  KW 


DO  43  ROW  =  1,PQ 
A(^ROW,l)  =  SI 


AC^ROW 

43  CONTINUE 


(ROW,l)  *  KW(1,1) 


WRITE  OUT  MEAN,  COVARIANCE,  AND  FILTER  COEFFICIENTS  TO  THE 
USER  INPUT  FILE. 

OPEN  (yNIT»2,FILE®FILTER( I), STATUS* 'NEW' ,CARRIAGECONTROL='LIST' , 
*  FORM* ' FORMATTED ' ) 
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c 

50 
500 

51 

52 
C 

C 

53 

54 
19 


WRITE 

FORMAT 


fl 


DO  51  J  a  1  PQ 

FORMAT  (F10.5) 

CX0SE(UN1T=2) 

WRITE (*,53) 
format'-  ■ 


WRITE 
FORMAT V. 
CONTINUE 
STOP 
END 


*  ^n^Tf9,^?f^^ICIENTS 


M2,$) 
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it  ic 

*  PROGRAM  SGMT2  * 

*  * 

*  PURPOSE  To  segment  a  single  image  with  two  textures  given  the  * 

*  image,  the  image  dimensions,  the  filter  dimensions,  * 

*  and  two  sets  ox  filter  parameters  to  be  used.  * 

*  REQUIRED  IMSL  ROUTINES  * 

*  NONE  * 

*  * 

*  IMPLEMENTED  BY  LTJG  TIMUR  KUPELI  Sep  1986  * 

it  it 

’kitit’k'kititititit’kititicit^ifkititit'k^ititit^ititititititititititic^iticitititicicit^^ititicicic^icicit^ititit'kitititicitit'kit 


INTEGER  IN,N,IM,P1,P,Q1,Q,TMAX,T,PQ,R0W,C0L,I,J,JJ,JJJ,L,LL, 
LLL , LLLL , COUNT , LI , HI , L J , HJ 


KW(2),MEAN(2) ,AA(4,1 ,2) ,TEMP,SUM1 ,SUM2 ,TEXTUR(128, 128, 2) , 
LN(2),ERR0R(128,128,2),PML1,PML2,PML(128,128), 
AREA,KS,A(2,2) 


CHARACTER*50 


IMAGE , FNAME , MLTEST , MAPTEST 


BYTE  BINPUT(128,128),ML(128,128),MAP(128,128),MLI(128,128) 

IN  =  128 
IM  =  128 
PI  =  4 
Q1  *  4 
TMAX  =  2 

GET  THE  INPUT  PARAMETERS  OF  THE  PROGRAM 
WRITE(*,2)  IN 

FORMAT ('  ENTER  THE  NUMBER  OF  ROWS  IN  IMAGE. LIMIT  OP' ,13 ,  ’  ; ' 

READ(*,3)  N 
FORMAT(I3) 

IF((N.LT.l)  .OR.  (N.GT.IN))  GOTO  1 
WRITE(*,5)  IM 

FORMAT('  ENTER  THE  NUMBER  OF  COLUMNS  IN  IMAGE. LIMIT  OF ' , 13 , ' : ' , 
READ(*,3)  M 

IF((M.LT.l)  .OR.  (M.GT.IM))  GOTO  4 
WRITE(*,7)  PI 

FORMATC  ENTER  THE  NUMBER  OF  ROWS  IN  FILTER. LIMIT  OF' ,13, ' : ' 
READ(*,3)  P 

IF((P.LT.2)  .OR.  (P.GT.Pl))  GOTO  6 
WRITE(*,9)  01 

FORMATC'  ENTER  THE  NUMBER  OF  COLUMNS  IN  FILTER. LIMIT  OF' ,13, '  i ' ,$) 
READ(*,3)  Q 

IF((Q.LT.2)  .OR.  (Q.GT.Ql))  GOTO  8 
WRITE  (*  11)  TMAX 

FORMaH*  ENTER  NUMBER  OF  TEXTURES  TO  PROCESS  .LIMIT  OF' ,13,  '  i '  ,$) 
REiU?(\3)  T  „  _  _  _ 


READ(*,3)  T 

IF((T.LT.2)  .OR.  (T.GT.TMAX))  GO  TO  10 

GET  THE  SINGLE-CHANNEL  IMAGE 

PQ  «  P  *  Q 
WRITE (*, 20) 

FORMATC  ENTER  FILENAME  OF  IMAGE  ' ,$) 
READ (^,21)  IMAGE 
FORMAT(A50) 
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OPEN (UNIT=1,FILE=IMAGE,STATUS=' OLD' ,ACCESS= ' DIRECT ' ) 

DO  22  ROW  =  1  .  N 

READ(l'ROW)  (BINPUT(ROW,COL),COL  =  1  ,  M  ) 

CONTINUE 

CLOSE (  UNIT  =  1  ) 

GET  THE  FILTER  COEFFICIENTS .MEANS ,  AND  COVARIANCES 

DO  23  I  =  1  ,  T 

WRITE (*  24)  I 

FORMATC  ENTER  THE'fILENAME  OF  FILTER  PARAMETERS  SET  NUMBER' ,12, 
READ(*,21)  FNAME 

0PEN(UNIT=2 , FILE=FNAME , STATUS® ' OLD ' ) 

READ(2,240)  KW(I) 

FORMAT(F10.5) 


FORMAT (FIO. 5) 
READ(2,240)  MEAN(I) 

DO  242  KK=  1,PQ 
READ(2,240)  AACKK,!,!] 


CONTINUE 
CLOSE (UNIT  =  2) 

WRITE(*,241)  KW(I),MEAN(I),AA(1,1,I),AA(2,1,I),AA(3,1,I),AA(4,1,I) 
FORMAT(6F10.5) 

CONTINUE 

READ (*,220)  MLTEST 
READ (*,220)  MAPTEST 
FORMAT (A80) 

CONVERT  A  2-D  BYTE, INPUT  IMAGE  DATA  IN  THE  RANGE  OF  -128  TO  127  THAT 
REPRESENTS  APPROPRIATE  INTENSITY  LEVELS  IN  THE  RANGE  OF  0  TO  255 

DO  30  ROW  ®  1  ,  N 

DO  31  COL  =  1  ,  M 


TEMP  »  BINPUT(ROW,COL) 
IF(TEMP.LT.O.O)  TEMP  *  TEMP+256 
TEXTUR(ROW,COL,l)  ®  TEMP  -  MEAN( 
TEXTUR(ROW,COL,2)  =  TEMP  -  MEAN( 


CONTINUE 


CONTINUE 


CALCULATOIN  OF  ERROR  ESTIMATE  FOR  TWO  TEXTURES 

DO  40  I  =  1  ,  T 
J  =  1 

DO  41  JJ  =  1  ,  P 

DO  42  JJJ  =  1  ,0 

A(JJ,JJJ)  =  ia(J,l,I) 

J  =  J  +  1 
CONTINUE 
CONTINUE 

DO  43  L  =  1  ,  N 
DO  44  LL  s  1  ,  M 

ERROR(L,LL,I)  =  0.0 
DO  45  LLL  =  1  ,  P 
J  =  L  -  LLL 
DO  46  LLLL  =1,0 

JJ  =  LL  -  LLLL 

IF((J.GT.O)  .AND.  (JJ.GT.O))  THEN 
ERROR(L,LL,I)  =ERROR(L,LL,I)+  A(LLL,LLLL)*TEXTUR(J, JJ,I) 
END  IF 

CONTINUE 


non 


CONTINUE 

CONTINUE 

CONTINUE 

LN(I)  *  ALOG(KW(I)) 

CONTINUE 

CALCULATION  OF  MAXIMUM  LIKELIHOOD  IMAGE  SEGMENTATION 

OPEN ( UNIT=3 , FILE=mTEST , STATUS* ' NEW '  , ACCESS* 'DIRECT ' 
RECL*(IM/4) ,MAXREC*IN) 

L  *  1 
LL  *  2 

DO  50  I  =  1  ,  N 
DO  51  J  =  1  ,  M 
PMLl  *  0.0 
PML2  =0.0 

PMLl  =  (  ERRORa,J,L)**2)/KW(L))  +  LN(L) 

IF{Pra,l  .GT.  PML2)  THEN 
ML(I,J)=  -1 
END  IF 

CONTINUE^''''"’  ' 


CLOSE (UNIT*3) 

MAXIMUM  A  POSTERIORI  IMAGE  SEGMENTATION 

,,CCESS='DIEECT'  , 

KS  *  100 

DO  600  II  *  1  ,  5 
DO  60  I  *  1  ,N  . 

DO  61  J  =  1  ,  M 
SUMl  =  0.0  ■ 

SUM2  *  0.0 
LI  *  I  -  3 
HI  *  I  +  3 
LJ  *  J  -  3 
HJ  *  J  +  3 

IF(LI.EQ.O)  LI  =  1 
IF(HI.GT.N)  HI  =  N 


IFfHI.GT.N 

IFfLJ.EQ.O 

IF(HJ,GT.M 


HI  =  N 
LJ  *  1 
HJ  =  M 


AREA  =  (HI  -  LI  +  1)  *  (HJ  -  LJ  +1) 

DO  62  ROW  =  LI  ,  HI 
DO  63  COL  =  LJ  ,  HJ 

IF((ROW.NE.I)  .OR.  (COL.NE.J))  THEN 

END  IF  ’  ”^I(R0W,C0L) 

CONTINUE 

CONTINUE 

MAP(I,J)*0 
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60 

600 


70 


ML1(1,J)  =  MAP(1,J) 
CONTINUE  ' 

CONTINUE 
CONTINUE 

DO  70  I  =  1  ,  N 

CONTINUE 

STOP 

END 
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